You made a small arithmetic error during your substitution, which when corrected should give you the correct result. Using the same substitution $u = 1 - 2xh + h^2 \implies \frac{\mathrm{d}u}{-2h} = \mathrm{d}x$
\begin{align*}
\int_{-1}^1 \frac{1}{1-2xh + h^2} \mathrm{d}x &= \frac{1}{-2h} \int_{(1+h)^2}^{(1-h)^2} \frac{\mathrm{d}u}{u} \\
&= -\frac{1}{2h}\log\frac{(1-h)^2}{(1+h)^2} \\
&= \frac{1}{h} (\log(1+h) - \log(1-h)) \\
&= \frac{1}{h} \left(\sum_{n=1}^\infty \frac{h^n}{n} - \frac{(-h)^n}{n}\right) \\
&= \frac{1}{h} \sum_{n\ \mathrm{odd}} \frac{2h^n}{n} \\
&= \frac{2}{h} \sum_{l = 0}^\infty \frac{h^{2l+1}}{2l+1} \\
&= \sum_{l = 0}^\infty \frac{2h^{2l}}{2l+1}
\end{align*}
Equating this to the right hand side, we have
$$\sum_{l = 0}^\infty \frac{2h^{2l}}{2l+1} = \sum_{l=0}^\infty h^{2l}\int_{-1}^1 P_l^2 \mathrm{d}x$$
Equating coefficients gives the desired result:
$$ \int_{-1}^1 P_n^2 \mathrm{d}x = \frac{2}{2n+1}$$
We know,from generating function of Legendre polynomials,
$(1-2xz+z^2)^\frac{-1}{2}=\sum z^nP_n(x)$
Squaring both sides and writing the sum by separating the cross terms and "self-quadrature" terms,we get
$(1-2xz+z^2)^{-1}=\sum z^{2n}P_n^2(x)+2\sum z^{m+n}P_m(x)P_n(x)$
Integrating both sides between -1 and +1,we have,
$\int_{-1}^1\sum z^{2n}P_n^2(x)dx+\int_{-1}^1 2\sum z^{m+n}P_m(x)P_n(x)dx=\int_{-1}^{1}(1-2xz+z^2)^{-1}dx$
Now,
$\int_{-1}^{1}P_n(x)P_m(x)dx=0$ due to orthogonality of Legendre polynomials and the 2nd integral on LHS of the above equation is wrt x,the terms containing m and n can be taken out,and therefore the 2nd integral vanishes.
$\Rightarrow \int_{-1}^1\sum z^{2n}P_n^2(x)dx+0=\int_{-1}^{1}(1-2xz+z^2)^{-1}dx$
The integral on RHS,i.e.,
$
{\displaystyle\int}\dfrac{1}{-2zx+z^2+1}\,\mathrm{d}x $ can be solved as--
Substitute $u=-2zx+z^2+1$
$\dfrac{\mathrm{d}u}{\mathrm{d}x} = -2z$
$\Rightarrow \mathrm{d}x=-\dfrac{1}{2z}\,\mathrm{d}u$
The integral reduces to
$=-\class{steps-node}{\cssId{steps-node-1}{\dfrac{1}{2z}}}{\displaystyle\int}\dfrac{1}{u}\,\mathrm{d}u$
which is a standard integral.
Therefore the integral becomes
$=-\dfrac{\ln\left(u\right)}{2z}$
$=-\dfrac{\ln\left(-2zx+z^2+1\right)}{2z}$
$\Rightarrow \sum z^{2n}\int_{-1}^1P_n^2(x)dx=\frac{-1}{2z}[\ln(1-2xz+z^2)]_{-1}^1$
RHS=
$\frac{-1}{2z}\ln\frac{1-2z+z^2}{1+2z+z^2}=\frac{-1}{2z}\ln(\frac{1-z}{1+z})^2=\frac{1}{z}[\ln(1+z)-\ln(1-z)]$
Using Maclaurin Series of $\ln(1+z)$ and $\ln(1-z)$ and expanding,
RHS=
$\frac{1}{z}[(z-\frac{z^2}{2}+\frac{z^3}{3}-\frac{z^4}{4}+...)-(-z-\frac{z^2}{2}-\frac{z^3}{3}-\frac{z^4}{4}-...)]$
Observe that the terms containing even powers of x cancel and the odd terms gets doubled(due to presence of them in both the series indicated above).Taking the 1/z term inside the square brackets,
$\Rightarrow \int_{-1}^1\sum z^{2n}P_n^2(x)dx=2[1+\frac{z^2}{3}+\frac{z^4}{5}+..\frac{z^{2n}}{2n+1}.]$
Finally,equating the coefficients of $z^{2n}$ on both sides we have,
$\int_{-1}^{1} (P_n(x))^2 dx = \frac{2}{2n+1}$,which is the required relation.
Best Answer
You recurrence rule seems to be wrong. What Wiki says is link: $$ P'_{l+1}(x)=(l+1)P_l(x)+xP'_l(x). $$ Then the computation becomes straightforward $$ \int_0^1 P_l(x) = \frac{1}{l+1}\int_0^1 (l+1)P_l(x) dx \\ = \frac{1}{l+1}\int_0^1 [P'_{l+1}(x)-xP'_l(x)]dx \\ = \frac{1}{l+1}\int_0^1 P'_{l+1}(x)dx - \frac{1}{l+1}\int_0^1 xP'_l(x)dx \\ = \frac{1}{l+1} [P_{l+1}(x)]_0^1 - \frac{1}{l+1}\left([xP_l]_0^1 - \int_0^1 P_l(x)dx\right) \\ = \frac{1}{l+1}[1-P_{l+1}(0)] - \frac{1}{l+1}\left(1 - \int_0^1 P_l(x)dx\right) \\ = \frac{1}{l+1}[1-P_{l+1}(0)] - \frac{1}{l+1} + \frac{1}{l+1} \int_0^1 P_l(x)dx \\ $$ We used integration by parts. Noting $I=\int_0^1 P_l(x)dx$ and $A=1-P_{l+1}(0)$ we need to solve the equation $$ (l+1)I= A-1+I \\ lI= A-1 \\ I = (A-1)/l\\ = -P_{l+1}(0)/l $$ There is a formula for value at zero: link