Well, the best way I can think of to introduce a real physical problem where Laguerre ànd Legendre (even of the associated form!) differential equations play a role, is the (exact, that is, without taking spin into consideration) solution of the combined angular momentum / magnetic moment / energy eigenvalue-problem for the hydrogen atom.
I cannot tell the whole story here, but I can sketch the process.
When solving the Schrodinger equation involved, and transforming it to spherical coordinates, all applicable quantum numbers roll into your lap and on route you will encounter the Laguerre and associated Legendre differential equations.
Now for the approach of solving these equations. Let's have a look at the Laguerre equation first.
$$xy'' + (1+\alpha-x)y' + ny = 0$$
If you try $y = \sum_{j=0}^{\infty}a_jx^j$, you will find that
$$\frac{a_{j+1}}{a_j} = -\frac{n-j}{(j+1)(j+1+\alpha)}$$
So, when moving from one coefficient to the next, the sign changes, it is divided by $j+1$, multiplied by $n-j$ and divided by $j+1+\alpha$. When j reaches n, the coefficients vanish.
That means that we will not be far off if we try
$$a_j = \frac{(-1)^j}{j!}\binom{n+\alpha}{n-j}$$
In fact, this is exactly what we need.
Now, which function will this be? The factor $\frac{(-1)^j}{j!}$ points at $e^{-x}$.
The factor $\binom{n+\alpha}{n-j}$ points at the $j^{th}$ derivative of $x^{n+\alpha}$, and we saw above that the series ends for j=n. But when we try $[e^{-x}x^{n+\alpha}]^{(n)}$, firstly we have to compensate for all the $e^{-x}$ factors by multiplying the result by $e^x$ , and secondly we would get terms for the running index j that have a factor $x^{n + \alpha - (n-j)} = x^{j + \alpha}$ and we would like to have $x^j$. So we have to compensate at the end here too, by multiplying by $x^{-\alpha}$. This reasoning leads to trying
$$y = e^xx^{-\alpha}[e^{-x}x^{n+\alpha}]^{(n)}$$
Where $f^{(n)}$ stands for the $n^{th}$ derivative of $f$.
And funnily enough, it's spot on again.
Over to the Legendre equation, first the ordinary one.
$$[(1-x^2)y']' + n(n+1)y = 0$$
Trying $y = \sum_{j=0}^{\infty}a_jx^j$ again, we find that
$$\frac{a_{j+2}}{a_j} = \frac{j(j+1) - n(n+1)}{(j+1)(j+2)}$$
As can be seen in one of my other posts (About the Legendre differential equation) there are solutions that are polynomials ($P_n$) and solutions that are represented by an unending series expansion ($Q_n$).
The $P_n$ are either odd-powered or even powered polynomials, depending on $n$ which is also the degree of $P_n$.
The $Q_n$ contain a factor $ln(\frac{1+x}{1-x})$, see (About the Legendre differential equation) again.
Finally, the associated Legendre equation.
$$[(1-x^2)y']' + \{l(l+1) - \frac{m^2}{1-x^2}\}y = 0$$
It is natural to try $y = (1-x^2)^{\alpha}P_l$. It is easily checked that we get the desired term $\frac{m^2}{1-x^2}y$ if we take $\alpha = m/2$. Further analysis shows that if we wish to arrive at all the correct terms, we have to use $P_l^{(m)}$ instead of $P_l$. Note that I still use the notation $f^{(m)}$ as denoting the $m^{th}$ derivative of $f$. The notation I use for the solutions of the associated Legendre equation is
$$P_{l,m} = (1-x^2)^{m/2}P_l^{(m)}$$
I hope this helps :-).
You have to multiply $f_n$ by $(1-x^2)^2$ to get an orthogonal family of polynomials. In general, if $\{P_n\}_{n\ge 1}$ is the system of orthogonal polynomials with respect to the weight $W(x)$ then $k-$ the derivatives form an orthogonal system with respect to the weight $W^k(x)a(x)$ where $b(x)$ comes from the differential equation for orthogonal polynomials:
$$a(x)y''(x)+b(x)y'+\lambda_n y=0.$$
In this particular case, Legandre polynomials arise as solutions of the differential equation
$$(1-x^2)y''(x)-2xy'+n(n+1)y=0.$$
and $W=1$ which makes $4- $th derivatives to be orthogonal with respect to the weight $(1-x^2)^4.$ In other words $(1-x^2)^2f_n$ are going to be orthogonal with respect to the weight $1.$
One way to prove the statement above is to use the given differential equation together with integration by parts.
Best Answer
The generating function for Legendre polynomials is \begin{eqnarray*} \sum_{n=0}^{\infty} P_n(z) h^n =\frac{1}{\sqrt{1-2hz+h^2}}. \end{eqnarray*} Differentiate this w.r.t $z$ \begin{eqnarray*} \sum_{n=0}^{\infty} P^{'}_n(z) h^n =\frac{(-1/2)(-2h)}{(1-2hz+h^2)^{3/2}}. \end{eqnarray*} Set $z=1$ \begin{eqnarray*} \sum_{n=0}^{\infty} \color{red}{P^{'}_n(1)} h^n =\frac{h}{(1-h)^{3}} =\sum_{n=0}^{\infty} \color{red}{\binom{n+1}{2}} h^n. \end{eqnarray*}