sequences and series - On twisted Euler sums

An interesting investigation started here and it showed that
$$ \sum_{k\geq 1}\left(\zeta(m)-H_{k}^{(m)}\right)^2 $$
has a closed form in terms of values of the Riemann $\zeta$ function for any integer $m\geq 2$.
I was starting to study the cubic analogue $ \sum_{k\geq 1}\left(\zeta(m)-H_{k}^{(m)}\right)^3 $ and I managed to prove through summation by parts that
$$ \sum_{k\geq 1}\frac{\left(H_k^{(2)}\right)^2}{k^2} =\frac{1}{3}\zeta(2)^3-\frac{2}{3}\zeta(6)+\zeta(3)^2 $$
where the LHS, according to Flajolet and Salvy's notation, is $S_{22,2}$. An explicit evaluation of $\sum_{k\geq 1}\left(\zeta(2)-H_{k}^{(2)}\right)^3 $ is completed by the computation of
$$\boxed{ S_{12,2} = \sum_{k\geq 1}\frac{H_k H_k^{(2)}}{k^2} }$$
which has an odd weight, hence it is not covered by Thm 4.2 of Flajolet and Salvy. On the other hand they suggest that by the kernel $(\psi(-s)+\gamma)^4$ the previous series and the cubic Euler sum $\sum_{n\geq 1}\frac{H_n^3}{(n+1)^2}=\frac{15}{2}\zeta(5)+\zeta(2)\,\zeta(3)$ are strictly related.

Question: can you help me completing this sketch, in order to get an explicit value for $S_{12,2}$ and for $\sum_{k\geq 1}\left(\zeta(2)-H_{k}^{(2)}\right)^3 $? Alternative techniques to summation by parts and residues are equally welcome.

Update: I have just realized this is solved by Mike Spivey's answer to Zaid's question here. On the other hand, Mike Spivey's approach is extremely lengthy, and I would be happy to see a more efficient derivation of $S_{12,2}=\zeta(2)\,\zeta(3)+\zeta(5)$.

Answer

I am skipping the derivation of Euler Sums $\displaystyle S(1,1;3) = \sum\limits_{n=1}^{\infty} \frac{H_n^2}{n^3}$ and $\displaystyle S(2;3) = \sum\limits_{n=1}^{\infty} \frac{H_n^{(2)}}{n^3}$ (for a derivation see solutions to $\textbf{problem 5406}$ proposed by C.I. Valean, Romania in SSMA).

Now, consider the partial fraction decomposition, \begin{align*}\sum\limits_{k=1}^{n-1}\left(\frac{1}{k(n-k)}\right)^2 = \frac{2}{n^2}\left(H_n^{(2)} + \frac{2H_n}{n}-\frac{3}{n^2}\right) \qquad \cdots (\star)\end{align*}

Multiplying both sides of $(\star)$ with $H_n$ and summing over $n \ge 1$ (and making the change of varible $m = n+k$):

\begin{align}
\sum\limits_{n=1}^{\infty}\frac{H_n}{n^2}\left(H_n^{(2)}+2\frac{H_n}{n} - \frac{3}{n^2}\right) &= \sum\limits_{n=1}^{\infty}\sum\limits_{k=1}^{n-1} \frac{H_n}{k^2(n-k)^2} \\&= \sum\limits_{m,k=1}^{\infty} \frac{H_{m+k}}{m^2k^2} \tag{0}\\&= \sum\limits_{m,k,j=1}^{\infty} \frac{mj+kj}{m^2k^2j^2(m+k+j)} \tag{1} \\&= 2\sum\limits_{m,k,j=1}^{\infty} \frac{jk}{m^2k^2j^2(m+k+j)} \tag{2} \\&= 2\sum\limits_{j,k=1}^{\infty} \frac{1}{jk}\sum\limits_{m=1}^{\infty} \frac{1}{m^2(m+j+k)} \tag{3} \\&= 2\sum\limits_{k,j=1}^{\infty} \frac{1}{kj(j+k)}\left(\zeta(2) - \frac{H_{k+j}}{k+j}\right) \tag{4} \\&= 4\sum\limits_{n=2}^{\infty}\frac{H_{n-1}}{n^2}\left(\zeta(2) - \frac{H_n}{n}\right)\\&= 4\zeta(2)\sum\limits_{n=1}^{\infty}\frac{H_{n-1}}{n^2} + 4\sum\limits_{n=1}^{\infty} \frac{H_n}{n^4}-4\sum\limits_{n=1}^{\infty} \frac{H_n^2}{n^3}\end{align}

Where, in steps $(0)$ and $(3)$ we used the identity, $\displaystyle \frac{H_q}{q} = \sum\limits_{\ell = 1}^{\infty} \frac{1}{\ell (\ell + q)}$. In step $(1)$ we used the symmetry w.r.t. the variables, in step $(2)$ interchanged order of summation and in step $(4)$ made the change of variables $n = j+k$.

Thus, $$ \sum\limits_{n=1}^{\infty} \frac{H_nH_n^{(2)}}{n^2} = 4\zeta(2)\zeta(3) + 7\sum\limits_{n=1}^{\infty} \frac{H_n}{n^4} - 6\sum\limits_{n=1}^{\infty} \frac{H_n^2}{n^3}$$