Proof orthogonality of Legendre polynomials using Sturm Liouiville theory only I want to prove that Legendre Polynomials corresponding to different eigenvalues are orthogonal, using (if possible) only Sturm Liouville Theory. I had a look here but they seem to prove something harder and with different techniques.
I will highlight my questions, so it is clearer where I have doubts.
$1$ GENERAL THEORY
$1.1$ General Sturm-Liouville Problem
A regular SL problem on $[x_0,x_1]$ is $$\frac{d}{dx}\left(p(x)\frac{dy}{dx}\right)+q(x)y+\lambda(x)y=0$$ with $p(x)$, $p'(x)$, $q(x)$ and $w(x)$ continuous on $[x_0,x_1]$ and with $p(x)>0$ and $w(x)>0$ on $[x_0,x_1]$.
$1.2$ General Sturm-Liouville Linear Operator
We have a linear operator $\mathcal L\colon y\mapsto \mathcal L[y]:=-\frac{1}{w(x)}\left(\frac{d}{dx}\left(p(x)\frac{dy}{dx}\right)+q(x)y\right)$ with the property that $\mathcal Ly=\lambda y$. Hence the SL equation is an eigenvalue problem.
$1.3$ Linear Operator is Self-Adjoint
After defining $$\langle f,g\rangle_w=\int_{x_0}^{x_1}f(x)g(x)w(x)\,dx$$ we say $\mathcal L$ is self-adjoint if for all pairs of functions $y_a$ and $y_b$ we have $\langle \mathcal L[y_a], y_b\rangle_w=\langle y_a, \mathcal L[y_b]\rangle_w$. By writing this in form of an integral, we obtain that it works, provided suitable boundary conditions (zero-Dirichlet, zero-Neumann,  radiation, singular points, periodic).
$1.4$ Properties of Solutions and Eigenvalues
We then then prove various things among which "Eigenfunctions corresponding to distinct eigenvalues are orthogonal".
$2$ LEGENDRE CASE STUDY
Legendre Differential Equation
Is given by $$\frac{d}{dx}\left((1-x^2)\frac{dy}{dx}\right)+n(n+1)y=0$$
hence setting $p(x)=1-x^2$, $q(x)=0$ amd $w(x)=1$ we have a SL equation, with $\lambda=n(n+1)$.
Legendre Equation as Sturm-Liouville Problem on $[-1,1]$
Now $w(x)=1>0$ however we also require $p(x)=1-x^2>0$. This is true for $x\in [-1,1]$. Notice also that $p(x)$, $p'(x)$, $q(x)$ and $w(x)$ are all continuous on this interval. Hence we have a singular Sturm Liouville problem on $[-1,1]$.
Legendre Equation Linear Operator
We have found a linear operator $\mathcal L\colon y \mapsto \mathcal L[y]:=-\frac{d}{dx}\left((1-x^2)\frac{dy}{dx}\right)$. This enables us to write $\mathcal Ly=n(n+1)y$ or $\mathcal Ly = \lambda y$ where $\lambda=n(n+1)$.

Can we just state that the Legendre polynomials (which are solutions to the Legendre equation) are orthogonal if they correspond to distinct eigenvalues (e.g.. $\lambda_n=n(n+1)$ and $\lambda_m=m(m+1)$ with $m\neq n$)?

Or do we have to prove it? If so, why do we have to prove it and how?
I mean it seems to me that we have the theory for the general Sturm Liouville problem. However in my textbook and in various pdfs found online, they prove everything saying that it is not sufficient to say it. How come?
 A: The Legendre equation is a singular equation on $(-1,1)$ because $p(x)=1-x^2$ vanishes at $x=\pm 1$. Without endpoint conditions, this operator is not well-defined as a selfadjoint Sturm-Liouville operator. The conditions required at $x=\pm 1$ are not obvious. There are singular types of Sturm-Liouville problems where endpoint conditions cannot be imposed, but that's not the case here. There are non-standard conditions at $x=\pm 1$ that lead to fully selfadjoint operators, with corresponding eigenfunctions that are not the Legendre polynomials, and for which the eigenvalues are not $n(n+1)$, as they are for the Legendre polynomials.
Without specifying the endpoint conditions at $x=\pm 1$, every real number is an eigenvalue, in the same way that regular problems do not have discrete eigenvalues without imposing endpoint conditions. As it turns out, the standard Legendre operator $Lf=-\frac{d}{dx}\left((1-x^2)\frac{df}{dx}\right)$ can be formulated by requiring that $f$ be twice continuously differentiable on $(-1,1)$, that $f$ be in $L^2$, that $Lf \in L^2(-1,1)$, and that $f$ remain bounded near $\pm 1$. That's enough to force a discrete set of eigenvalues--namely $\lambda=n(n+1)$--and to give orthogonality of the eigenfunctions. Without such conditions, though, this is not a well-posed Sturm-Liouville problem. For example, $Lf=0\cdot f$ has two independent $L^2$ solutions $f_1 = 1$, $f_2 = \ln\frac{1+x}{1-x}$. And $f_2$ is not orthogonal to the polynomial solution $p(x)=x$ of $Lp=2p$, even though $x$ is orthogonal to $f_1=1$. So, it is absolutely true that something must be said about endpoint conditions, in other to obtain a selfadjoint problem where eigenfunctions corresponding to different eigenvalues are automatically orthogonal.
The Legendre equation is in the "limit-circle" case at $x=\pm 1$, which means that all classical solutions of $Lf=\lambda f$ are in $L^2(-1,1)$, regardless of $\lambda\in\mathbb{C}$. And that means that endpoint conditions are needed at each endpoint in order to obtain a selfadjoint problem, which is quite analogous to the case for regular Sturm-Liouville problems on a finite interval. The endpoint conditions, however, are non-trivial to describe. For example, if $f$ is twice continuously differentiable in $(-1,1)$ with $Lf \in L^2$, then there exist unique numbers $A_{\pm}(f)$ and $B_{\pm}(f)$ such that
$$
          f(x) = A_{\pm}(f)\frac{1}{2}\ln\left(\frac{1+x}{1-x}\right)+B_{\pm}(f)+o\left(\sqrt{1-x^{2}}\ln\left(\frac{1+x}{1-x}\right)\right), \;\;\; x\approx  \pm 1.
$$
The numbers $A_{\pm}(f)$, $B_{\pm}(f)$ are computed as limits:
\begin{align}
          A_{\pm}(f) & =  \lim_{x\rightarrow \pm 1}\left[(1-x^{2})f'(x)\right], \\
          B_{\pm}(f) & = \lim_{x\rightarrow\pm 1}\left[f(x)-A_{\pm}(f)\frac{1}{2}\ln\left(\frac{1+x}{1-x}\right)\right].
\end{align}
Valid endpoint conditions at $x=\pm 1$ that lead to selfadjoint problems are analogous to ordinary endpoint conditions, and have the general form
$$
          \cos\alpha A_{-}(f) + \sin\alpha B_{-}(f) = 0 \\
          \cos\beta A_{+}(f) + \sin\beta B_{+}(f) = 0.     \tag{$\dagger$}
$$
The conditions giving rise to the ordinary Legendre polynomials as eigenfunctions are:
$$
       A_{-}(f) = 0,\;\; A_{+}(f) = 0.   \tag{$\dagger\dagger$}
$$
Imposing conditions $(\dagger\dagger)$ forces $f$ to have limits at $\pm 1$.
But there are selfadjoint problems for $L$ that do not result in the Legendre polynomials as solutions; a two parameter family of such operators $L_{\alpha,\beta}$ are described by $(\dagger)$. It turns out that the above classical conditions $(\dagger\dagger)$ are equivalent to imposing a condition of boundedness on $f$ near $\pm 1$, which is the classical way that Physicists and Engineers generally pose the problem, not realizing that this is an actual endpoint condition where some endpoint quantity is $0$. The functionals $A_{\pm},B_{\pm}$ as described above form a basis for the continuous linear functionals on the graph of $L$ that vanish on $\mathcal{C}_{c}^{\infty}(-1,1)$, which is how abstract boundary functionals are defined for singular Sturm-Liouville operators.
A: This is perhaps easier to understand from linear algebra. Consider the eigenvalue problem with symmetric matrix $A$ (e.g. $A = A^T$). Let $(\lambda_1,u)$ and $(\lambda_2,v)$ be eigenvalue/vector pairs. Consider
$$ v^T A u - u^T A v = (\lambda_1 - \lambda_2) \langle u, v \rangle. $$
Finally note that $A = A^T$ implies that the LHS is 0 and so either $\lambda_1 = \lambda_2$ or the eigenvectors are orthogonal. The SL operator is a symmetric linear operator and so shares this property.
