Integral identity with square of Jacobi polynomial This has stumped me for a while: I have a function $\zeta_k^S(x)$ that can be expressed using Jacobi polynomials $P_k^{(\alpha,\beta)}(x)$:
$\zeta_k^S(x):=\displaystyle\sum_{j=0}^k\frac{(-1)^jx^{2j+2}{k \choose j}(4S-k-1)!}{(j+2)!(1+x^2)^{2S}(4S-k-j-1)!}=\frac{x^2k!}{\Gamma(3+k)(1+x^2)^{2S}}P_k^{(2,-4S-2)}(1+2x^2)$
where $S\in\{3,\frac{7}{2},4,\frac{9}{2},\ldots\}$ and $k\in\{0,1,\ldots,2S-1\}$ and  using 
$\frac{\Gamma(m-n)}{\Gamma(n)}=\Bigg\{\begin{array}{cc} (-1)^m\frac{n!}{(n-m)!} & : m\leq n\\ 0 & : m>n \end{array}\ \ \ \ \ \ \ \ \ $ for $n\in\{-1,-2,-3,\ldots\}\ \ $.
I want to show the following identity, which I believe to be true: 
$\mathcal{I}_k :=8\pi S\displaystyle\int_0^{\infty}\frac{\zeta_k^S(x)\cdot\zeta_k^S(x)\cdot x\ dx}{(1+x^2)^2}=\frac{4\pi S}{(4S-2k-1)(4S-k)(4S-k+1)(k+1)(k+2)}\ \ $.
But I'm at a loss for how. Here's what I have tried:
$\bullet$ Doing the integral explicitly. I get the following sum:
$\mathcal{I}_k=\displaystyle\sum_{j,m=0}^k\frac{4\pi S(-1)^{j+m}{k \choose j}{k \choose m}(4S-k-1)!^2(j+m+2)!(4S-j-m-2)!}{(j+2)!(m+2)!(4S-k-j-1)!(4S-k-m-1)!(4S+1)!}\ \ $.
I've tried to simplify this and also to use proof by induction, but I can't seem to relate $\mathcal{I}_{k+1}$ to $\mathcal{I}_k$ in a useful way.
$\bullet$ Using Rodriguez' formula $P_k^{(\alpha,\beta)}(x)=\frac{(-1)^k}{2^kk!}(1-x)^{-\alpha}(1+x)^{-\beta}\frac{d^k}{dx^k}\Big((1-x)^{\alpha+k}(1+x)^{\beta+k}\Big)$ and the expression for the derivative of the Jacobi polynomial, together with partial integration. I make some headway, but fail in the end.
$\bullet$ Again using proof by induction, but this time through the following recurrence relation 
$2n(n+\alpha+\beta)(2n+\alpha+\beta-2)P_n^{(\alpha,\beta)}=(2n+\alpha+\beta-1)\Big((2n+\alpha+\beta)(2n+\alpha+\beta-2)x-\alpha^2-\beta^2\Big)P_{n-1}^{(\alpha,\beta)}-2(n-\alpha-1)(n+\beta-1)(2n+\alpha+\beta)P_{n-2}^{(\alpha,\beta)}$
This almost works, but I haven't found a way to deal with the extra factor of $x$ in front of $P_{n-1}^{(\alpha,\beta)}$, and this messes up my relation between $\mathcal{I}_k$, $\mathcal{I}_{k-1}$ and $\mathcal{I}_{k-2}$.
$\bullet$ Looking at recurrence relations between hypergeometric functions, in which $\zeta_k^S(x)$ can also be expressed, without finding any useful ones.
$\bullet$ Using the expression in this paper relating products of Jacobi functions to linear combinations. However it ends up being too hairy for me. 
Any ideas?
 A: Let me first rewrite the statement in a bit more understandable form. 

Under the restrictions you have on $S$ and $k$, you want to show that 
$$\int_0^{\infty}\frac{\left[P_k^{(2,-4S-2)}(1+2x^2)\right]^2x^5dx}{(1+x^2)^{2+4S}}=\frac{(k+1)(k+2)}{2(4S-2k-1)(4S-k)(4S-k+1)}.\tag{1}$$

To compute the integral on the left, rewrite it as
\begin{align}
\int_0^{\infty}\frac{\left[P_k^{(2,-4S-2)}(1+2x^2)\right]^2x^5dx}{(1+x^2)^{2+4S}}=\frac{1}{2}\int_0^{\infty}\frac{\left[P_k^{(2,-4S-2)}(1+2y)\right]^2y^2dy}{(1+y)^{2+4S}}=\\
=2^{4S-2}\int_1^{\infty}(1-t)^2(1+t)^{-4S-2}\left[P_k^{(2,-4S-2)}(t)\right]^2dt,\tag{2}
\end{align}
where we used the change of variables $y=x^2$ to obtain the 1st equality and $t=1+2y$ for the 2nd one. Note that (2) looks very much like $\|P_k^{(2,-4S-2)}\|^2$ with respect to the usual orthogonality measure $w(t)=(1-t)^2(1+t)^{-4S-2}$ for Jacobi polynomials $P_k^{(2,-4S-2)}$, except that we integrate from $1$ to $\infty$ instead of from $-1$ to $1$. Actually, the latter integral would diverge because of values of $S$ but this hints nevertheless that (2)  should be a "good" quantity.
Indeed, let us replace $w(t)P_k^{(2,-4S-2)}(t)$ in (2) using Rodrigues formula for Jacobi polynomials:
$$w(t)P_k^{(2,-4S-2)}(t)=\frac{1}{(-2)^k k!}\frac{d^k}{dt^k}\left(1-t\right)^{2+k}\left(1+t\right)^{-4S-2+k}$$
Integrate $k$ times by parts (it is rather clear that all the boundary terms will vanish) so that all derivatives act on the second polynomial in the product:
\begin{align}&\int_1^{\infty}(1-t)^2(1+t)^{-4S-2}\left[P_k^{(2,-4S-2)}(t)\right]^2dt=\\
&=\frac{1}{2^k k!}
\int_1^{\infty}\left(1-t\right)^{2+k}\left(1+t\right)^{-4S-2+k}\frac{d^k}{dt^k}\left(P_k^{(2,-4S-2)}(t)\right)dt=\\&=\frac{A_k}{2^k}\int_1^{\infty}\left(1-t\right)^{2+k}\left(1+t\right)^{-4S-2+k}dt,\tag{3}
\end{align}
where $A_k$ denotes the leading coefficient in $P_k^{(2,-4S-2)}(t)=A_kt^k+\ldots$ This coefficient is explicitly given by
$$A_k=\frac{(4S-k-1)!}{(-2)^k k!\,(4S-2k-1)!}.\tag{4}$$
In addition, the integral in (3) can also be explicitly calculated (e.g. setting $u=\frac{t-1}{t+1}$ one obtains a beta function integral):
$$\int_1^{\infty}\left(1-t\right)^{2+k}\left(1+t\right)^{-4S-2+k}dt=(-1)^k2^{-4S+1+2k}
\frac{(4S-2k-2)!\,(k+2)!}{(4S-k+1)!}.\tag{5}$$
Now combining (2), (3), (4) and (5), after careful elementary simplification we obtain (1). $\blacksquare$
