<s>[INST] Challenging Sum: compute $\sum_{n=1}^\infty\frac{H_n}{2n+1}\left(\zeta(3)-H_n^{(3)}\right)$ 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.
 [/INST] A second solution in large steps by Cornel Ioan Valean
Let's start with the following useful identity which is easily derived by using recurrence relations and simple rearrangements, manipulations with sums, that is

Let $n$ be a non-negative integer number. Then, we have
  $$\int_0^1 x^{2n}\frac{\log(1+x)}{1+x}\textrm{d}x$$
$$=\frac{1}{2}H_{2n}^2-2\log(2) H_{2n}+\frac{1}{2}H_{2n}^{(2)}-\frac{1}{4}H_n^2-\frac{1}{4}H_n^{(2)}+\log (2)H_n+\frac{1}{2} \log ^2(2)-\sum_{k=1}^{n-1}\frac{H_k}{2 k+1},$$
  where $H_n^{(m)}=1+\frac{1}{2^m}+\cdots+\frac{1}{n^m}$ represents the $n$th generalized harmonic number of order $m$.

By multiplying both sides of the identity above by $1/n^3$ and considering the summation from $n=1$ to $\infty$, we get
$$\sum_{n=1}^{\infty} \frac{1}{n^3}\sum_{k=1}^{n-1}\frac{H_{k}}{2 k+1}=\sum_{k=1}^{\infty} \sum_{n=k+1}^{\infty}\frac{1}{n^3}\frac{H_{k}}{2 k+1}=\underbrace{\sum_{k=1}^{\infty}\frac{H_{k}}{2 k+1}\left(\zeta(3)-H_k^{(3)}\right)}_{\text{The desired series}}$$
$$=\frac{1}{2}\sum_{n=1}^{\infty}\frac{H_{2n}^2}{n^3}-2\log(2) \sum_{n=1}^{\infty}\frac{H_{2n}}{n^3}+\frac{1}{2}\sum_{n=1}^{\infty}\frac{H_{2n}^{(2)}}{n^3}-\frac{1}{4}\sum_{n=1}^{\infty}\frac{H_n^2}{n^3}-\frac{1}{4}\sum_{n=1}^{\infty} \frac{H_n^{(2)}}{n^3}$$
$$+\log (2)\sum_{n=1}^{\infty} \frac{H_n}{n^3}+\frac{1}{2}\log ^2(2)\sum_{n=1}^{\infty}\frac{1}{n^3}-\int_0^1 \frac{\log(1+x)}{1+x}\operatorname{Li}_3(x^2)\textrm{d}x,$$
where we see all the series in the right-hand side are easily reducible to known series which may also be found in the book (Almost) Impossible Integrals, Sums, and Series. 
On the other hand, with simple integration by parts, we obtain
$$\int_0^1 \frac{\log(1+x)}{1+x}\operatorname{Li}_3(x^2)\textrm{d}x$$
$$=\frac{1}{2}\log^2(2)\zeta(3)-2\int_0^1 \frac{\log^2(1+x)\operatorname{Li}_2(x)}{x}\textrm{d}x-2\int_0^1 \frac{\log^2(1+x)\operatorname{Li}_2(-x)}{x}\textrm{d}x,$$
where the last integrals may be found calculated in the paper The calculation of a harmonic series with a weight $5$ structure, involving the product of harmonic numbers, $H_n H_{2n}^{(2)}$.
A note: The sister of the result above (easy to obtain by recurrence relations and very useful), 

$$\int_0^1 x^{2n-1} \frac{\log(1+x)}{1+x}\textrm{d}x$$
$$=2\log(2) H_{2n}-\log(2)H_n+\frac{1}{4}H_n^2+\frac{1}{4}H_n^{(2)}-\frac{1}{2}H_{2n}^2-\frac{1}{2} H_{2n}^{(2)}+\frac{H_{2n}}{2n}-\frac{H_n}{2n}
$$
$$
-\frac{1}{2}\log^2(2)+\sum_{k=1}^{n-1}\frac{H_k}{2 k+1}.
$$

</s>