How to prove
$$\sum_{n=1}^\infty\frac{H_n}{2n+1}\left(\zeta(3)-H_n^{(3)}\right)=\frac74\zeta(2)\zeta(3)-\frac{279}{16}\zeta(5)+\frac43\ln^3(2)\zeta(2)-7\ln^2(2)\zeta(3)\\+\frac{53}4\ln(2)\zeta(4)-\frac2{15}\ln^5(2)+16\operatorname{Li}_5\left(\frac12\right)$$
where $H_n^{(q)}=\sum_{k=1}^n\frac{1}{k^q}$ is the generalized harmonic number, $\operatorname{Li}_a(x)=\sum_{k=1}^\infty\frac{x^k}{k^a}$ is the polylogarithmic function and $\zeta$ is the Riemann zeta function.
This problem was proposed by Cornel and no solution has been submitted yet. I managed to convert it to a double integral but it seems tough to crack. Here is what I did:
Using the integral representation of the polygamma function:
$$\int_0^1\frac{x^n\ln^a(x)}{1-x}dx=-\psi^{(a)}(n+1)=(-1)^a a!\left(\zeta(a+1)-H_n^{(a+1)}\right)$$
With $a=2$ we have
$$\zeta(3)-H_n^{(3)}=\frac12\int_0^1\frac{x^n\ln^2(x)}{1-x}dx\overset{x=y^2}{=}4\int_0^1\frac{y^{2n+1}\ln^2(y)}{1-y^2}dy$$
multiply both sides by $\frac{H_n}{2n+1}$ then sum up we get
$$\sum_{n=1}^\infty\frac{H_n}{2n+1}\left(\zeta(3)-H_n^{(3)}\right)=4\int_0^1\frac{\ln^2(y)}{1-y^2}\left(\sum_{n=1}^\infty\frac{y^{2n+1}H_n}{2n+1}\right)dy$$
we have
$$\sum_{n=1}^\infty \frac{y^{2n+1}H_n}{2n+1}=-\int_0^y\frac{\ln(1-x^2)}{1-x^2}dx$$
which follows from integrating $\sum_{n=1}^\infty x^{2n}H_n=-\frac{\ln(1-x^2)}{1-x^2}$ from $x=0$ to $x=y$.
so
$$\sum_{n=1}^\infty\frac{H_n}{2n+1}\left(\zeta(3)-H_n^{(3)}\right)=-4\int_0^1\int_0^y\frac{\ln^2(y)\ln(1-x^2)}{(1-y^2)(1-x^2)}dxdy$$
$$=-4\int_0^1\frac{\ln(1-x^2)}{1-x^2}\left(\int_x^1\frac{\ln^2(y)}{1-y^2}dy\right)dx$$
For the inner integral, Mathematica gives
$$\int_x^1\frac{\ln^2(y)}{1-y^2}dy\\=\operatorname{Li}_3(-x)-\operatorname{Li}_3(x)-\ln(x)\operatorname{Li}_2(-x)+\ln(x)\operatorname{Li}_2(x)-\ln^2(x)\tanh^{-1}(x)+\frac74\zeta(3)$$
and the integral turned out very complicated. So any good idea how to approach the harmonic series or the integral?
Thank you.