This post was prompted by this question at Mathematics Stack Exchange. This proof of the orthogonality of the Legendre polynomials is from *Special Functions and Their Applications* by N. N. Lebedev, a book that I highly recommend.

We begin with Legendre’s differential equation

\begin{equation}

[(1-x^{2})P^{\prime}_{n}(x)]^{\prime} +n(n+1)P_{n}(x) = 0,\quad n \in \mathbb{Z}_{0}^{+}

\label{eq:lp1}

\tag{1}

\end{equation}

The first step is to multiple equation \eqref{eq:lp1} by \(P_{m}(x)\) and subtract it from equation \eqref{eq:lp1} written for \(m\) and multiplied by \(P_{n}(x)\).

\begin{equation}

[(1-x^{2})P^{\prime}_{m}(x)]^{\prime}P_{n}(x) \,-\, [(1-x^{2})P^{\prime}_{n}(x)]^{\prime}P_{m}(x) + [m(m+1)-n(n+1)]P_{m}(x)P_{n}(x) = 0

\end{equation}

Rearrangement yields

\begin{equation}

\{(1-x^{2})[P^{\prime}_{m}(x)P_{n}(x)-P^{\prime}_{n}(x)P_{m}(x)]\}^{\prime} + (m-n)(m+n+1)P_{m}(x)P_{n}(x) = 0

\label{eq:lp2}

\tag{2}

\end{equation}

Integrating equation \eqref{eq:lp2} from -1 to 1, the first term goes to 0 and we are left with

\begin{equation}

(m-n)(m+n+1) \int\limits_{-1}^{1} P_{m}(x)P_{n}(x) dx = 0

\end{equation}

or

\begin{equation}

\int\limits_{-1}^{1} P_{m}(x)P_{n}(x) dx = 0, \quad m \ne n

\end{equation}