Precise conditions on Sturm-Liouville Theorems In Sturm-Liouville (SL) theory (https://en.wikipedia.org/wiki/Sturm-Liouville_theory), there are three fundamental theorems concerning the solutions of the SL differential equation,
$ \frac{\mathrm{d}}{\mathrm{d}x}\left[p(x)\frac{\mathrm{d}y(x)}{\mathrm{d}x}\right]+q(x)y(x)=\lambda y(x)$.
They are (SL Theorems):


*

*The set of eigenvalues $\lambda=\{\lambda_1, \lambda_2,\ldots\}$ are all real, countable and distinct.

*The set of eigenfunctions $y(x)=\{y_1(x),y_2(x),\ldots\}$, forms a orthogonal system in some interval $(a,b)$ of the real line, so that they satisfy $\int_a^b y_m(x)y_n(x)\mathrm{d}x=K_n\delta_{m,n}$, where $K_n$ are non-null constants. 

*This set of eigenfunctions forms a basis for the vector space of square integrable functions. 


Now, it is very known that a second-order differential equation of the form
$f(x)y''(x)+g(x)y'(x)+h(x)y(x)=\lambda y(x)$
can be put into the SL form after we multiply it by some integrating function $w(x)$ that satisfies the first-order differential equation
$\frac{\mathrm{d}}{\mathrm{d}x}\left[w(x)f(x)\right] = w(x)g(x)$.
In this case, the solutions of this differential equation will obey the SL theorems, if we agree that the solutions are orthogonal with respect to the weight function $w(x)$, i.e., they satisfy $\int_a^b y_m(x)y_n(x)w(x)\mathrm{d}x=K_n\delta_{m,n}$.
My question is: what are the precise conditions (on the functions $p,q,f,g,h,w$, on the interval of orthogonality etc.) for the SL Theorems above to hold?
I ask this motivated by the following related problem: it is very known that some second-order differential equations admit an infinite sequence of orthogonal polynomials on the real line -- these are called the Classical Orthogonal Polynomials (COP) and comprehends the Jacobi, Hermite and Laguerre polynomials). There is a theorem due to Bochner (see for instance the book of T. Chihara, "An Introduction to Orthogonal Polynomials", p. 150) that these three sequences of COP are the only ones satisfying a second-order differential equation that are orthogonal on the real line (up to linear transformations). However, in Bochner proof appears another infinite sequence of polynomials -- called nowadays as Bessel polynomials -- that although are not orthogonal on the real line, they are orthogonal on the complex unit circle. These Bessel polynomials satisfy the differential equation (see Krall, H.L., and Frink, O. "A new class of orthogonal polynomials: The Bessel polynomials." Transactions of the American Mathematical Society 65.1 (1949): 100-115.): 
$x^2y''(x) + 2(x+1)y'(x)=n(n+1)y(x)$.
It can be verified that $w(x)=\exp(-2/x)$ is the integrating function for this differential equation (hence the weight function), so it can be put into the SL form. However, I can't see why the SL Theorems do not hold for it (they seem to not hold, since, for instance, the Bessel polynomial of degree 2, $B_2(x)=3x^2+3x+1$, has imaginary roots, while the roots of any orthogonal polynomial sequence on the real line are all real). 
I will appreciate any comment about these questions. 
 A: The SL theorem you cite can be applied to the operator
$$
        Lf = \frac{1}{w}\left[-\frac{d}{dx}\left(p\frac{df}{dx}\right)+q\right]
$$
on $[a,b]$ if $w$ is positive and continuous on $[a,b]$, if $p$ is positive and continuously differentiable on $[a,b]$, and if $q$ is continuous on $[a,b]$. The theorem is then applied to the eigenvalue problem with a fixed set of endpoint conditions:
$$
             Lf=\lambda f,\\
        \cos\alpha f(a)+\sin\alpha f'(a)=0\\
        \cos\beta f(b)+\sin\beta f'(b)=0.
$$
(Here $\alpha,\beta$ are real and fixed, and up to you to chose.) In this case you get a complete orthogonal set of eigenfunctions, and the eigenspaces are one-dimensional. The space you have to work in is $L^2_w[a,b]$, where
$$
       \langle f,g\rangle = \int_{a}^{b}f(x)g(x)w(x)dx.
$$
There are various ways to extend this. You can extend this theory to cases where $p,w$ vanish at $a$ or at $b$, or where the interval is infinite or semi-infinite, provided all classical solutions of $Lf=0$ are in $L^2_w[a,b]$, which is a condition where both endpoints are regular or in the so-called "limit-circle case" (this is standard terminology in the field.) However, in the extended cases, the endpoint conditions must be altered in a non-obvious way that I won't go into. Basically, there are some asymptotic expansions near $a,b$ that take the place of endpoint values. I posted an example of the Legendre operator
$$
       Lf = \frac{d}{dx}\left((1-x^2)\frac{df}{dx}\right),\;\; -1 < x < 1.
$$
This equation is in the limit circle case at both endpoints because all solutions of $Lf=0$ are in $L^2(-1,1)$ (the weight function here is $1$.) It turns out that boundedness near both endpoints--which is the classical condition--is a valid asymptotic condition for this equation, and you end up with a nice selfadjoint problem in that case. However, the associated Legendre equations are not in the limit-circle case; so no endpoint conditions can be imposed, and you still end up with discrete eigenfunctions, even though the theory I mentioned above does not apply to the associated Legendre equations.
Selfadjoint Restrictions of Legendre Operator $-\frac{d}{dx}(1-x^{2})\frac{d}{dx}$
