How to prove
$$\sum_{n=1}^\infty\frac{H_{n}^2}{(2n+1)^3}=\frac{31}{8}\zeta(5)-\frac{45}{8}\ln2\zeta(4)+\frac72\ln^22\zeta(3)-\frac78\zeta(2)\zeta(3)$$
where $H_n$ is the harmonic number and $\zeta$ is the Riemann zeta function.
Here is my approach and would like to see different methods if possible.
Using the identity
$$\frac{\ln^2(1-x)}{1-x}=\sum_{n=1}^\infty (H_n^2-H_n^{(2)})x^n$$
replace $x$ with $x^2$, then multiply both sides by $\frac12\ln^2x$ and integrate from $x=0$ to $1$ we get
$$\sum_{n=1}^\infty\frac{H_n^2-H_n^{(2)}}{(2n+1)^3}=\frac12\int_0^1\frac{\ln^2x\ln^2(1-x^2)}{1-x^2}\ dx\\=\frac1{16}\int_0^1\frac{\ln^2x\ln^2(1-x)}{\sqrt{x}(1-x)}\ dx=\frac1{16}\left.\frac{\partial^4}{\partial a^2\partial b^2}\text{B}(a,b)\right|_{a\mapsto 1/2\\b\mapsto0^{+}}$$
with help of Mathematica we have
$$\left.\frac{\partial^4}{\partial a^2\partial b^2}\text{B}(a,b)\right|_{a\mapsto 1/2\\b\mapsto0^{+}}=248\zeta(5)-90\ln2\zeta(4)+56\ln^22\zeta(3)-112\zeta(2)\zeta(3)$$
Also, from here we have
$$\sum_{n=1}^\infty\frac{H_{n}^{(2)}}{(2n+1)^3}=\frac{49}{8}\zeta(2)\zeta(3)-\frac{93}{8}\zeta(5)$$
Combining these two result, we get the desired answer.
Note: You can find here details about the derivative of beta function.