This is from here
$$\int^1_{-1}f_n(x)P_n(x)dx = 2(-1)^n\frac{a_n}{2^n}\int^1_0(x^2-1)^ndx=2(-1)^n\frac{a_n}{2^n}.I_n$$........(6)
I don't understand as in shouldnt it be like this, $$\int^1_{-1}f_n(x)P_n(x)dx = (-1)^n\frac{a_n}{2^n}\int^1_{-1}(x^2-1)^ndx=0$$ as they should cancel out even if the integral is non-zero.
Edited: Lastly, how does $\int^1_{-1}f_n(x)P_n(x)dx$ shows orthogonality?
Here is a much simpler proof from Special Functions and Their Applications by N. N. Lebedev.
We begin with Legendre's differential equation \begin{equation} [(1-x^{2})P^{\prime}_{n}(x)]^{\prime} +n(n+1)P_{n}(x) = 0,\quad n \in \mathbb{Z}_{0}^{+} \label{eq:lp1} \tag{1} \end{equation}
The first step is to multiple equation \eqref{eq:lp1} by $P_{m}(x)$ and subtract it from equation \eqref{eq:lp1} written for $m$ and multiplied by $P_{n}(x)$.
\begin{equation} [(1-x^{2})P^{\prime}_{m}(x)]^{\prime}P_{n}(x) \,-\, [(1-x^{2})P^{\prime}_{n}(x)]^{\prime}P_{m}(x) + [m(m+1)-n(n+1)]P_{m}(x)P_{n}(x) = 0 \end{equation}
Rearrangement yields \begin{equation} \{(1-x^{2})[P^{\prime}_{m}(x)P_{n}(x)-P^{\prime}_{n}(x)P_{m}(x)]\}^{\prime} + (m-n)(m+n+1)P_{m}(x)P_{n}(x) = 0 \label{eq:lp2} \tag{2} \end{equation}
Integrating equation \eqref{eq:lp2} from -1 to 1, the first term goes to 0 and we are left with \begin{equation} (m-n)(m+n+1) \int\limits_{-1}^{1} P_{m}(x)P_{n}(x) dx = 0 \end{equation} or \begin{equation} \int\limits_{-1}^{1} P_{m}(x)P_{n}(x) dx = 0, \quad m \ne n \end{equation}