My attempt:
$P_m$ and $P_n$ satisfy the Legendre equation, so
(1)
\begin{equation} \label{1} (1-x^2)P_m''-2xP_m'+m(m+1)P_m=0 \end{equation}
(2)
\begin{equation} (1-x^2)P_n''-2xP_n'+n(n+1)P_n=0 \end{equation}
Multiplying (1) by $P_n'$ and (2) by $P_m'$ and adding and simplifying gives me:
$\frac{d}{dx} [(1-x^2) P_n' P_m']=2xP_n'P_m'+m(m+1)P_mP_n'+n(n+1)P_nP_m'$
I see that integrating from $-1$ to $1$ twice will give me desired LHS, but I can't prove RHS is $0$. Please help!
You have the right idea (using the definition by ODE) but are applying it in the wrong place. Let's instead integrate by parts: $$\int_{-1}^1 (1-x^2)P_m'P_n'dx = (1-x^2)P_m'P_n \Bigr |_{-1}^1 - \int_{-1}^1 [(1-x^2)P_m']'P_ndx = -\int_{-1}^1 [(1-x^2)P_m']'P_ndx $$ Then we have $$ [(1-x^2)P_m']' = (1-x^2)P_m'' - 2xP_m' = -m(m+1)P_m$$ $$\implies -\int_{-1}^1 [(1-x^2)P_m']'P_ndx = \int_{-1}^1 m(m+1)P_mP_ndx = 0$$ by orthogonality.