I am trying to prove that $\frac{1}{n!\space2^n}\frac{d^n}{dx^n}\{(x^2-1)^n\}=P_n(x)$, where $P_n(x)$ is the Legendre Polynomial of order n.
I've been told that the proof uses complex analysis, of which I know nothing, isn't there a proof with elementary methods? (If there isn't, I'm still interested in the other one).
If $n$ is integer \begin{eqnarray*} \frac{d^n}{d x^n} (x^2-1)^n &=& \frac{d^n}{d x^n} \left [ \sum_{k=0}^n (-1)^k \frac{n!}{k!(n-k)!} x^{2n-2k} \right ] \\ &=& \sum_{k=0}^n (-1)^k \frac{n!}{k! (n-k)!} \frac{(2n-2k)!}{(n-2k)!} x^{n-2k}. \end{eqnarray*} The sum above does not go up to $k=n$, since after $k=[n/2]$, the derivatives are 0 then we write
\begin{eqnarray*} \frac{d^n}{d x^n} (x^2-1)^n &=& \sum_{k=0}^{[n/2]} (-1)^k \frac{n!}{k! (n-k)!} \frac{(2n-2k)!}{(n-2k)!} x^{n-2k}. \end{eqnarray*} It follows from the infinite series truncated to the Legendre polynmial that
\begin{eqnarray*} P_n(x) = \frac{1}{2^n n!} \frac{d^n}{dx^n} (x^2-1)^n. \end{eqnarray*}
The approach followed here is in reverse order. We started with Rodriguez's formula and showed that it corresponds to a Legendre polynomial. A more intuitive approach is to start at the polynomials
\begin{eqnarray*} y(x)= (1-x^2)^n. \end{eqnarray*} and take derivates, and verifty that the derivatives taken $n$ times will get you to the Legendre differential equation. That is, we have that
\begin{eqnarray*} y' = -2 n x (1-x^2)^{n-1} \end{eqnarray*} which we can write as
\begin{eqnarray} (1-x^2) y' + 2n x y = 0. \label{tole} \end{eqnarray} and starts looking a bit like a Legendre differential equation.
We want to differentiate this equation $k$ times and use the Leibniz rule. That is, if we call $u=1-x^2$,
\begin{eqnarray*} \frac{ d^k}{dx^k} [u y'] = \sum_{j=0}^{k} \binom{k}{j} u^{(j)} y^{(k-j+1)} \end{eqnarray*} Given that $u$ is a second order polynomial only three terms of this sum will survive. That is
\begin{eqnarray*} \frac{ d^k}{dx^k} [u y'] &=& u y^{(k+1)} + k u' y^{(k)} + k(k-1) u^{(2)} y^{(k-1)} \\ &=& (1-x^2)y^{(k+1)} - 2 k x y^{(k)} -2 \frac{k(k-1)}{2} y^{(k-1)} = 0 \end{eqnarray*} Likewise we use the Leibniz rule for the product $2nxy$ where only two terms will survive. That is
\begin{eqnarray*} \frac{ d^k}{dx^k} [2 n x y] &=& 2 n x y^{(k)} + 2 n k y^{(k-1)}, \end{eqnarray*} we combine the two results above to find
\begin{eqnarray*} (1-x^2)y^{(k+1)} - 2 k x y^{(k)} - k(k-1) y^{(k-1)} + 2 n x y^{(k)} + 2 n k y^{(k-1)} = 0 \end{eqnarray*} At this point we observe that if $k=n+1$, we find
\begin{eqnarray*} (1-x^2)y^{(n+2)} - 2(n+1) x y^{(n+1)} - n(n+1) y^{(n)} + 2 n x y^{(n+1)} + 2 n (n+1) y^{(n)} = 0 \end{eqnarray*} which simplifies to
\begin{eqnarray*} (1-x^2) y^{(n+2)} - 2 x y^{(n+1)} + n(n+1) y^{(n)}=0. \end{eqnarray*}
and this is the Legendre differential equation with $y^n=P_n$. We then showed that
\begin{eqnarray*} \frac{d^n}{dx^n}(1-x^2)^n \end{eqnarray*} satisfies the Lagrange differential equation. The factor $1/(2^n n!)$ is included to make $P(1)=1$.
Check that the left side indeed defines an $n$th order polynomial.
Check that $\displaystyle \int_{-1}^{1}P_n(x)P_m(x)dx$ vanishes for $m\neq n$ (integration by parts).
Check the normalization condition $P_n(1)=1$ (Leibniz rule).
Added: As you almost correctly write in the comment below, the result of integration by parts (assuming that $m<n$ and transferring the derivatives from $P_n$ to $P_m$) can be written as $$\int_{-1}^1P_m(x)P_n(x)dx=\sum_{k=1}^{n}c_{mnk}\left[\frac{d^{m+k-1} (x^2-1)^m}{dx^{m+k-1}}\frac{d^{n-k}(x^2-1)^n}{dx^{n-k}}\right]_{-1}^{1},$$ where $c_{mnk}$ is some irrelevant constant. Consider the second factor in the square brackets. There you have a polynomial having $n$th order zeros at $x=\pm 1$ which we differentiate $n-k$ times. The result will therefore have $k$th order zeros at these points, which implies vanishing of the integral.
(The general formula of Legendre Polynomials is given by following equation: $$ P_k(x)=\sum_{m=0}^{\frac{k}{2}|\frac{k-1}{2}}{\frac{(-1)^m(2k-2m)!}{2^km!(k-m)!}}\frac{1}{(k-2m)!} x^{k-2m} $$
The Rodrigues' formula is:
$$ \frac{1}{2^kk!}\frac{d^k}{dx^k}[(x^2-1)^k] $$
The Binomial theorem is as follow: $$(x+y)^k=\sum_{i=0}^{k}\frac{k!}{i!(k-i)!}x^{k-i}y^{i}$$
Then $$(x^2-1)^k=\sum_{i=0}^{k}\frac{k!}{i!(n-i)!}{(x^{2})}^{k-i}(-1)^{i}$$
$$\frac{1}{2^kk!}\frac{d^k}{dx^k}[(x^2-1)^k]=\frac{1}{2^kk!}\frac{d^k}{dx^k}\big[\sum_{i=0}^{k}\frac{k!}{i!(n-i)!}{(x^{2})}^{k-i}(-1)^{i}\big]$$
So $$ =\frac{1}{2^kk!}\sum_{i=0}^{k}\frac{k!}{i!(n-i)!}\frac{d^k}{dx^k}(x)^{2k-2i}(-1)^{i} .....(1)$$
$$\frac{d^k}{dx^k}(x)^{2k-2i}=\frac{(2k-2i)!}{k-2i}x^{k-2i}......(2)$$
By compensate (2) into (1), we get:
$$ P_k(x)=\sum_{m=0}^{\frac{k}{2}|\frac{k-i}{2}}{\frac{(-1)^i(2k-2i)!}{2^ki!(k-i)!}}\frac{1}{(k-2i)!} x^{k-2i} $$
Change dummy variable (i) to (m), we get the same general formula of Legendre Polynomials :
$$ P_k(x)=\sum_{m=0}^{\frac{k}{2}|\frac{k-1}{2}}{\frac{(-1)^m(2k-2m)!}{2^km!(k-m)!}}\frac{1}{(k-2m)!} x^{k-2m} $$
I still think the other answers are a little too complicated. Our goal is to find an $n+1$-degree polynomial $P_{n+1}(x)$ on the interval $[a,b]$ that is orthogonal to any $n$-degree polynomial $Q_{\leq n}(x)$ on the interval $[a,b]$, i.e. $$\int_a^b P_{n+1}(x) Q_{\leq n}(x) dx = 0.$$ The LHS is an integral of a product, so hopefully one's first instinct would be to try integration by parts. Allow me to use the notation of the tabular method: https://kevinboone.me/parts.html (though I can't typeset the arrows here, even if it's possible with more technology) . We have 2 options: we can put $P_{n+1}$ in the derivative column, or the integral column. Let us look at putting it in the derivative column first: $$\begin{array}{r|c|c} \text{row }\# & \frac{d}{dx} & \int \\ \hline 0 & \color{blue}{P_{n+1}(x)} & Q_{\leq n}(x) \\ 1 & P_{n+1}'(x) & \color{blue}{\int Q_{\leq n}}\\ \vdots &\vdots & \vdots \\ n+1 & \color{red}{\text{constant}} & \ldots \\ n+2 & 0 & \color{red}{\ldots} \end{array}$$ The left column ("row $0$") starts with degree $n+1$, decreasing in degree by 1 each time, until on row $n+1$ it becomes constant (degree $0$). The last term in the integration by parts is then some [constant] times some high degree antiderivative of $Q_{\leq n}(x)$, which we have no control over in making $0$ upon evaluation $(\bullet)|^{x=b}_{x=a}$.
Let us try the other option, putting $P_{n+1}$ in the integral column: $$\begin{array}{r|c|c} \text{row }\# & \frac{d}{dx} & \int \\ \hline 0 & \color{blue}{Q_{\leq n}(x)} & P_{n+1}(x) \\ 1& Q_{\leq n}'(x) & \color{blue}{\int P_{n+1} }\\ \vdots & \vdots & \vdots \\ n& \color{red}{\text{constant}} & \ldots \\ n+1& 0 & \color{red}{(\int)^{n+1} P_{n+1}} \end{array}$$
In the scenario where $Q_{\leq n}$ has degree $n$ (on row $0$ of the table), then on row $n$, it will be degree $0$, namely a constant; on the right/integral column, on row $n+1$ we will have the $(n+1)$th antiderivative of $P_{n+1}(x)$. So, the result of this integral by parts table is (at most) $n$ terms: $$\int_a^b P_{n+1}(x) Q_{\leq n}(x) dx = \textstyle \Big[Q_{\leq n}(x) \cdot \int P_{n+1}(x)\Big]\bigg|_{x=a}^{x=b} - \ldots + (-1)^n \Big[\text{const} \cdot (\int)^{n+1} P_{n+1}(x)\Big]\bigg|_{x=a}^{x=b}$$ (The notation is a bit fiddly, but I hope you can see that the ideas are actually quite simple. If I was writing on a whiteboard in front of you, I could just draw slope-$(-1)$ ovals in the table representing each term, without having to use all this notation, but oh well.)
Again, recall that we have control over what $P_{n+1}$ is, and we have just found that all it needs to be is a degree $n+1$ polynomial s.t. its $1$st, $2$nd, ..., $(n+1)$th antiderivative all vanish at both $x=a$ and $x=b$.
This is still a little hard to think about, but it becomes much easier, if I tell you to think about the problem backward: start with a polynomial $A(x)$ of degree $2(n+1)$ (supposed/destined to be the $(n+1)$th antiderivative of $P_{n+1}$), so that its $0$th, $1$st, ..., $n$th derivatives vanish at $x=a$ and $x=b$ (and then taking $P_{n+1}(x):= (\frac d{dx})^{n+1} A(x)$, we're done!)
And indeed, hopefully one is familiar (from working with Taylor series, or orders of zeroes and poles in complex analysis) that the derivative of a function vanishing to order $M$ at a point $x=x_0$ vanishes to order $M-1$ at $x=x_0$. So, taking $A(x)=(x-a)^{n+1} (x-b)^{n+1}$ (indeed of degree $2(n+1)$), the $0$th, $1$st, ..., $n$th derivatives vanish at $x=a$ and $x=b$, as desired (and indeed the $(n+1)$th derivative is a polynomial of degree $n+1$, as desired!)
