Properties of solutions to the associated Legendre ODE The associated Legendre ODE is given by 
$$  \left( (1-x^2) f'(x) \right)' - \frac{m^2}{1-x^2} f(x) = \lambda f(x)$$
The eigenfunctions have certain properties that I would like to understand by looking NOT at the eigenfunctions and eigenvalues but only by referring to the ODE itself. So please pretend that you are not aware of the fact that you can construct the two things analytically.
Then, we have that the eigenfunctions are all $C^{\infty}$ if we exclude the Legendre-polynomials of the second kind, can we see this directly from the differential equation? Probably this follows for the Legendre polynomials somehow if we require boundedness of the solutions at the end of the interval to exclude the Legendre polynomials of the second kind. But how does this follow for the associated Legendre polynomials, i.e. then the ODE has an additional singular term that could cause additional trouble? Despite, all the solutions are proper well-behaved functions, so we should somehow see this out of the ODE itself.
Furthermore, if we increase $m$, then the ground-state solutions get higher. I mean, if we start with $m=0$, then the ground-state is $\lambda=0$, for $m=1$ we have $\lambda= 1(1+1)=2$ and for $m=2$ we have $\lambda = 2(2+1)=6$. Can we see this directly out of the differential equations that this is the case?  
 A: 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.
