Find this sum $S$ using Real-analysis methods only $$S = \sum_{k=1}^{\infty}\frac{2H_k}{(k+1)(k+2)^3}$$
I have tried a lot and failed, any help is appreciated. $H_k$ is the harmonic number.
Thanks (real method only please)
 A: We have:
$$\frac{1}{2}\log^2(1-x)=\sum_{k=1}^{+\infty}\frac{H_k}{k+1}x^{k+1}$$
but since $\int_{0}^{1} x^n\log^2 x\,dx = \frac{2}{(n+1)^3}$ it follows that:
$$\begin{eqnarray*} S &=& \frac{1}{2}\int_{0}^{1} \log^2(1-x)\log^2 x\,dx=\frac{1}{2}\left.\frac{\partial^4}{\partial^2 a\,\partial^2 b}B(a,b)\,\right|_{a,b=1}\end{eqnarray*}.$$
Since $B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}$ and $\Gamma'(z)=\psi(z)\Gamma(z)$ by differentiating multiple times we get a long expression depending only on the values of $\psi,\psi',\psi''$ and $\psi'''$ in $1$ and $2$:

$$ S = 12 - 4\,\zeta(2) -4\,\zeta(3)-\frac{1}{2}\,\zeta(4).$$


Addendum: An interesting approach is given by writing $\log(x)\log(1-x)$ as a linear combination of shifted Legendre polynomials. From the Rodrigues formula, exploiting integration by parts, we get:
$$\log x\log(1-x)=2-\frac{\pi^2}{6}-\sum_{j=1}^{+\infty}\frac{4j+1}{2j^2(2j+1)^2}P_{2j}(2x-1)$$
hence the Parseval's identity gives:
$$ S = \frac{1}{2}\left(\left(2-\frac{\pi^2}{6}\right)^2+\sum_{j=1}^{+\infty}\frac{(4j+1)}{4j^4(2j+1)^4}\right) $$
and the problem boils down to showing that:
$$\sum_{j=1}^{+\infty}\frac{(4j+1)}{4j^4(2j+1)^4}=20-4\zeta(2)-8\zeta(3)-\frac{7}{2}\zeta(4).\tag{1}$$
This problem seems to be feasible to creative telescoping techniques:
$$\frac{1}{(2n)^4}-\frac{1}{(2n+1)^4}=\frac{1+4n}{16n^4(2n+1)^4}+\frac{1+4n}{4n^3(2n+1)^3},$$
$$\frac{1}{(2n)^3}-\frac{1}{(2n+1)^3}=\frac{1}{8n^3(2n+1)^3}+\frac{3}{4n^2(2n+1)^2},$$
$$\frac{1}{(2n)^2}-\frac{1}{(2n+1)^2}=\frac{4n+1}{4n^2(2n+1)^2}=\frac{1}{4n^2}-\frac{1}{(2n+1)^2}.$$
