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 :-).
Best Answer
The equation must be considered in the context of the Hilbert space $L^{2}(-1,1)$ because the filter of requiring eigenfunctions to be in $L^{2}$ is what eliminates the non-regular solutions, and it is what determines the eigenvalues. For $m = 1,2,3,\cdots$, the operators $$ L_{m}f = -\frac{d}{dx}\left[(1-x^{2})\frac{d}{dx}f\right] + \frac{m^{2}}{1-x^{2}}f $$ are selfadjoint on the domain $\mathcal{D}(L_{m})$ consisting of all twice absolutely continuous functions $f$ for which $L_{m}f \in L^{2}$. No endpoint conditions are needed or are possible. So you can integrate by parts and be assured that all of the evaluation terms vanish in order to obtain $$ (L_m f,f) = \|\sqrt{1-x^{2}}f'\|^{2}+\|\frac{m}{\sqrt{1-x^{2}}}f\|^{2}. $$ It is automatically true that $f \in \mathcal{D}(L_{m})$ implies $$ \sqrt{1-x^{2}}f',\; \frac{1}{\sqrt{1-x^{2}}}f \in L^{2}(-1,1). $$ Hence, the product of these two expressions is in $L^{1}(-1,1)$, which gives $ff' \in L^{1}$ and thereby guarantees the existence of the following endpoint limits $$ \left.\int_{-1}^{1}ff'\,dx = \frac{f^{2}(x)}{2}\right|_{-1}^{1}. $$ These endpoint limits are also $0$ because there are non-zero boundary functionals; or you can appeal to the fact that $f/\sqrt{1-x^{2}}\in L^{2}$ in order to conclude that $f^{2}(\pm 1)$ cannot be non-zero. So quite a lot can be said just knowing that $f \in \mathcal{D}(L_m)$. And you can carry the analysis further by assuming $f \in \mathcal{D}(L_m)$ is also a solution of $L_m f = \lambda f$ for some $\lambda$.
Another important part of classical analysis for solving ODEs is the Method of Frobenius http://en.wikipedia.org/wiki/Frobenius_method . This classical method gives series solutions for equations with regular singular points, which nearly nearly all of the classical Sturm-Liouville eigenvalue problems have, at least in the finite plane. For example, $x=0$ is a regular singular point of $$ p(x)y''+q(x)y'+r(x)y = 0 $$ if $p$, $q$, $r$ have power series expansions around $x=0$, and the normalized equation $$ y''+\frac{q}{p}y+\frac{r}{p}y = 0 $$ has no worse than an order 1 pole for $q/p$ and an order 2 pole for $r/p$. Then you can get an approximation for the behavior of at least one solution by solving Euler's equation $x^{2}y''+axy'+by=0$ where $a$ and $b$ are the coefficients of the highest order singular terms of $q/p$ and $r/p$, respectively. This leads to at least one solution of the form $x^{m}\sum_{n=0}^{\infty}a_{n}x^{n}$ where $m$ is the solution of the indicial equation $m(m-1)+am+b=0$ with the largest real part. More generally, the substitution $y=x^{m}w$ leads to a new equation for $w$ which has a power series solutions at $x=0$. So that's the classical method that was invented about 140 years ago.
The power of the Method of Frobenius in this case is two-fold:
It motives a substitution $f=(1-x^{2})^{m/2}g$ which greatly simplifies the equation, and where it can be directly seen that differentiatng a solution $y=g$ for some $m$ and $\lambda$ gives a solution $y=g'$ for $m+1$ and the same $\lambda$.
The substitution leads to an equation which admits power series solutions. Because of the symmetry at $\pm 1$, the new equation admits entire solutions which happen to be polynomials for specific $\lambda$. A direct power series analysis shows that only the polynomial ones are acceptable solutions in $\mathcal{D}(L_m)$.
To carry out the Method of Frobenius, start with the eigenfunction equation: $$ (1-x^{2})f''-2xf'-\frac{m^{2}}{1-x^{2}}f+\lambda f = 0 \\ (x^{2}-1)f''+2xf'-\frac{m^{2}}{x^{2}-1}f-\lambda f = 0 \\ f''+\frac{2x}{(x-1)(x+1)}f'-\left[\frac{m^{2}}{(x-1)^{2}(x+1)^{2}}+\frac{\lambda}{(x-1)(x+1)}\right]f = 0. $$ (Note: I have negated your eigenvalue parameter because $L_{m}$ is a positive operator; so $L_{m}f=\lambda f$ leads to $\lambda > 0$ using the negative of your $\lambda$.) Only the highest order terms are initially considered in this method. For example, consider the equation near $x=1$: $$ f'' + \left[\frac{1}{x-1}+\cdots\right]f'+\left[-\frac{m^{2}}{4(x-1)^{2}}+\cdots\right]f = 0. $$ This determines a form of solution $f=(x-1)^{\alpha}g$ where $\alpha$ satisfies the indicial equation $$ \alpha(\alpha-1)+\alpha - \frac{m^{2}}{4} = 0 \\ %% \alpha^{2}-\frac{m^{2}}{4} = 0,\\ \alpha = \pm \frac{m}{2}. $$ Because the difference of these roots is an integer, only the one with the largest real part (i.e., $\alpha=m/2$) is guaranteed in general to lead to a solution of the form $$ f(x)= (1-x)^{m/2}\sum_{n=0}^{\infty}a_{n}(1-x)^{n}. $$ So classical considerations suggest a substitution of the form $$ f(x) = (1-x^{2})^{m/2}g(x). $$ This substitution leads to a simpler equation which is also sometimes called the Associated Legendre Equation: $$ (1-x^{2})g''-2x(m+1)g'-m(m+1)g + \lambda g = 0. $$ I believe that it was this form of the equation where it was discovered that differentiating the equation lead to another equation of the same form, but with a different $m$. For example, differentiate once and you get a new equation in $h=g'$: $$ (1-x^{2})h''-2x(m+2)h'-(m+1)(m+2)h + \lambda h = 0. $$ You'll notice that the new equation is the same as the original, but with $m$ replaced by $m+1$. So I think you can see how taking derivatives of solutions of the base Legendre equation where $m=0$, $$ (1-x^{2})g''-2xg'+\lambda g = 0 $$ leads to solutions of all of the higher order equations. All you have to do is multiply the derivatives by factors $(1-x^{2})^{m/2}$ in order to obtain full solutions of your original equation. Explicitly, if $P_{n}$ is the Legendre polynomial of order $n$, then $$ (1-x^{2})^{m/2}\frac{d^{m}}{dx^{m}}P_{n}(x) $$ is a solution of your equation.