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.
The following is an adapted translation of Wikipedia.fr, Suite_de_polynômes_orthogonaux#Relation_de_récurrence.
Let us prove that
$$p_{n+1}=(x+b_n)p_n-c_np_{n-1}$$
where
$$b_n=k_{n+1}-
k_n,\qquad c_n=\frac{h_n}{h_{n-1}},$$
$k_j,h_j$ being given by
$$p_j(x)=x^j+k_jx^{j-1}+\cdots,\quad h_j=\langle p_j,p_j \rangle.$$
(By convention, $c_0=p_{-1}=k_0=0.$)
The polynomial
$(x+b_n)p_n–p_{n+1}$ has a degree $<n$ hence may be written $\sum_{j=0}^{n-1}\mu_{n,j}p_j,$
with
$h_j\mu_{n,j}=\langle (x+b_n)p_n-p_{n+1},p_j\rangle=\langle xp_n,p_j\rangle$
(since for $j<n,$ $p_j$ is orthogonal to $p_n$ and to $p_{n+1}$).
Moreover,
$\langle xp_n,p_j\rangle=\langle p_n,xp_j\rangle.$
- For $j<n-1,$ this inner product is $0$ because $\deg(x p_j)<n.$
- For $j=n-1,$ it is equal to $h_n$ because (reasoning like previously) $x p_{n-1}–p_n$ has a degree $<n.$
We may conclude:
$(x+b_n)p_n-p_{n+1}= c_np_{n-1},$
where $c_n=\mu_{n,n-1}=\frac{h_n}{h_{n-1}}.$
Best Answer
Note that the left-hand term in your last equation is $0$, since $((1-x^2)^n)^{(n-1)} = 0$ for $x = -1, 1$. This follows because $(1-x^2)^n$ has degree-$n$ zeroes at $x = 1, -1$ and, in general, when a polynomial $p$ has a degree-$n$ zero at $x_0$ then $p(x_0), p'(x_0), \ldots, p^{(n-1)}(x_0) = 0$.
More generally, $((1-x^2)^n)^{(n-k)} = 0$ for $x = -1, 1$ for any $k > 0$. Repeatedly using this integration by parts trick - i.e., using induction - you can show that
$$\int_{-1}^1 Q_m(x)((1-x^2)^n)^{(n)} dx = (-1)^k \int_{-1}^1 Q_m^{(k)}(x)((1-x^2)^n)^{(n-k)} dx$$
for any $k$. Then when you reach $k = m+1$, you see that $Q_m^{(k)}(x) = 0$ since $Q_m$ is a polynomial of degree $m$, and the result follows.