Integral involving Legendre polynomial and $x^n$ I am trying to show that,
\begin{align}
I = \int_{-1}^1 x^nP_n(x)\,\mathrm{d}x = \frac{2^{n+1}n!n!}{(2n+1)!}
\end{align}
So far I have done the following. Rodrigues formula is as follows:
\begin{align*}
P_n(x) = \sum_{k=0}^N \frac{(-1)^k (2n-2k)!}{2^nk!(n-k)!(n-2k)!} x^{n-2k}
\end{align*}
where,
\begin{align*}
\begin{aligned}
N&=n/2, && \text{if} \quad n=\text{even} \\
N&=(n-1)/2, && \text{if} \quad n=\text{odd}
\end{aligned}
\end{align*}
Substitute Rodrigues formula,
\begin{align*}
I
&=
\int_{-1}^1 x^n \sum_{k=0}^N \frac{(-1)^k (2n-2k)!}{2^nk!(n-k)!(n-2k)!} x^{n-2k} \,\mathrm{d}x \\
&=
\sum_{k=0}^N \frac{(-1)^k (2n-2k)!}{2^nk!(n-k)!(n-2k)!} \int_{-1}^1 x^{2n-2k} \,\mathrm{d}x \\
&=
\sum_{k=0}^N \frac{(-1)^k (2n-2k)!}{2^nk!(n-k)!(n-2k)!} \left. \frac{x^{2n-2k+1}}{2n-2k+1}\right\rvert_{-1}^1 \\
&=
\sum_{k=0}^N \frac{(-1)^k (2n-2k)!}{2^nk!(n-k)!(n-2k)!} \frac{1 - (-1)^{2n-2k+1}}{2n-2k+1} \\
&=
\sum_{k=0}^N \frac{(-1)^k (2n-2k)!}{2^nk!(n-k)!(n-2k)!} \frac{1 + (-1)^{2n-2k}}{2n-2k+1}
\end{align*}
Since $2n-2k$ is even,
\begin{align*}
I
&=
\sum_{k=0}^N \frac{(-1)^k (2n-2k)!}{2^nk!(n-k)!(n-2k)!} \frac{2}{2n-2k+1} \\
&=
\sum_{k=0}^N \frac{2^{1-n}(-1)^k (2n-2k)!}{k!(n-k)!(n-2k)!(2n-2k+1)}
\end{align*}
Can someone give me a hint about how to proceed or do I need to prove with another way?

It is faster to exploit the generating function that comes from Rodrigues' formula or Bonnet's recursion formula:
$$\frac{1}{\sqrt{1-2xt+t^2}} = \sum_{n=0}^\infty P_n(x) t^n\tag{1} $$
By replacing $t$ with $xt$ we have:
$$ \int_{-1}^{1}\frac{dx}{\sqrt{1-2x^2 t+x^2 t^2}}=\sum_{n\geq 0}\left(\int_{-1}^{-1}x^n P_n(x)\,dx\right)t^n\tag{2} $$
so our integrals can be computed through the Taylor series of the LHS, regarded as a function of $t$.
The LHS of $(2)$ is:
$$ \frac{2\arcsin\!\sqrt{t(2-t)}}{\sqrt{t(2-t)}}=\frac{4\arcsin\sqrt{\frac{t}{2}}}{\sqrt{t(2-t)}} \tag{3} $$
and your initial claim turns out to be equivalent to a well-known identity.
Truth to be told, it is also very interesting sideways: some series related with the $\Gamma^2$ function can be easily computed by exploiting the mentioned property of Legendre polynomials.

Using Rodrigues' formula we have:
$$\int_{-1}^{1}x^n P_n(x)\,dx = \frac{1}{2^n n!}\int_{-1}^{1}x^n\cdot\frac{d^n}{dx^n}(x^2-1)^n\,dx \stackrel{RIBP}{=}\frac{1}{2^n}\int_{-1}^{1}(1-x^2)^n\,dx\tag{4}$$
where $\text{RIBP}$ stands for repeated integration by parts.
But the RHS of $(4)$ is a value of the Euler's beta function:
$$\frac{1}{2^n}\int_{-1}^{1}(1-x^2)^n\,dx=\frac{1}{2^n}\int_{0}^{1}z^{-1/2}(1-z)^n\,dz = \frac{\Gamma\left(\frac{1}{2}\right)\Gamma(n+1)}{2^n\,\Gamma\left(n+\frac{3}{2}\right)}\tag{5}$$
and we are equally done, since:
$$\frac{\Gamma\left(\frac{1}{2}\right)\Gamma(n+1)}{2^n\,\Gamma\left(n+\frac{3}{2}\right)} = \frac{2\cdot n!}{(2n+1)!!} = \frac{2^{n+1} n!^2}{(2n+1)!}.\tag{6}$$
