The generating function for Legendre polynomials is given by:
$$ g(x,t)=\frac{1}{\sqrt{1+t^2-2xt}}=\sum_{n\geq 0}P_n(x)\, t^n\tag{1} $$
hence integrating both sides over $(0,1)$ with respect to $x$,
$$\forall t\in(0,1),\qquad \frac{\sqrt{1+t^2}-1+t}{t}=\sum_{n\geq 0}\left(\int_{0}^{1}P_n(x)\,dx\right)t^n\tag{2} $$
and to compute our coefficients it is enough to compute the Taylor series of $\sqrt{1+t^2}$ around $t=0$. That proves formula $(49)$ linked by Omran Kouba and:
$$ f(x) = \frac{1}{2}+\sum_{m\geq 0}\frac{4m+3}{4m+4}\cdot\frac{(-1)^m}{4^m}\binom{2m}{m}\cdot P_{2m+1}(x).\tag{3} $$
Equality has to be intended as convergence in $L^2(-1,1)$. Gibbs phenomenon arises: it would be integersting to compute the magnitude of the overshoot near zero.
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
$$\sum_{n=0}^\infty \left(P'_n(x)\color{red}{t^n}-2xP'_n(x)t^{n+1}+P'_n(x)t^{n+2}\right)-\sum_{n=0}^\infty P_n(x)t^{n+1}=0$$ $$\color{red}{P_0'(x) + t(P_1'(x) -2xP_0'(x))~+} \sum_{\color{red}{n=2}}^\infty \left(P'_n(x)-2xP'_{n-1}(x) + P'_{n-2}(x)\right)t^n - \sum_{\color{red}{n=1}}^\infty P_{n-1}(x)t^n = 0$$ $$\color{red}{P_0'(x) + t(P_1'(x) -2xP_0'(x) - P_0(x))~+} \sum_{\color{red}{n=2}}^\infty \left(P'_n(x)-2xP'_{n-1}(x) + P'_{n-2}(x)-P_{n-1}\right)\color{red}{t^n} = 0$$
This is a power series in $t$. The only way a convergent power series can be $0$ for all $t$ in its radius of convergence is if all of its coefficients are $0$ (as can be proven by repeatedly differentiating with respect to $t$). Therefore $$P_0'(x) = 0\\P_1'(x) - 2xP_0'(x) - P_0(x) = 0$$ and for all $n \ge 2$, $$P'_n(x)-2xP'_{n-1}(x) + P'_{n-2}(x)-P_{n-1} = 0$$ All that remains is another shift of $n$ (which could have been avoided by being more careful earlier).