Compute $\sum_{n=1}^\infty\frac{H_n^2H_n^{(2)}}{n^3}$ How to prove

$$\sum_{n=1}^\infty\frac{H_n^2H_n^{(2)}}{n^3}=\frac{19}{2}\zeta(3)\zeta(4)-2\zeta(2)\zeta(5)-7\zeta(7)\ ?$$
  where $H_n^{(p)}=1+\frac1{2^p}+\cdots+\frac1{n^p}$ is the $n$th generalized harmonic number of order $p$.

This series is very advanced and can be found evaluated in the book (Almost) Impossible Integrals, Sums and Series page 300 using only series manipulations, but luckily I was able to evaluate it using only integration, some harmonic identities and results of easy Euler sums.
Can we prove the equality above in different methods besides series manipulation and the idea of my solution below? All approaches are highly appreciated. 
Solution is posted in the answer section. 
Thanks
 A: The series $\sum_{n=1}^\infty\frac{H_n^2H_n^{(2)}}{n^3}$ can be written as
$$\sum_{\substack{n_1\geq n_2\geq 1 \\ n_1\geq n_3\geq 1 \\ n_1\geq n_4\geq 1}}\frac{1}{n_1^3 n_2 n_3 n_4^2},$$
which can be recognized as a linear combination of multiple zeta values of weight $7$.
Multiple zeta values of weight $w$ are series of the form $$\zeta(s_1, \ldots, s_k) = \sum_{n_1 > n_2 > \cdots > n_k > 0} \ \frac{1}{n_1^{s_1} \cdots n_k^{s_k}},$$
such that $s_1,\dots,s_k$ are positive integers and $s_1>1$ such that $s_1+\dots+s_k=w$.
By breaking your sum into parts (depending whether $n_1>n_2>n_3>n_4$ or $n_1>n_2>n_3=n_4$ etc), your sum is equal to the following expression:
\begin{align*}
\sum_{n=1}^\infty\frac{H_n^2H_n^{(2)}}{n^3}&=2\zeta(3,2,1,1)+2\zeta(3,1,2,1)+2\zeta(3,1,1,2)+2\zeta(5,1,1)+2\zeta(4,2,1)+2\zeta(4,1,2)
\\&\quad
+\zeta(3,3,1)+2\zeta(3,2,2)+2\zeta(3,1,3)+2\zeta(6,1)+2\zeta(5,2)+2\zeta(4,3)
\\&\quad
+\zeta(3,4)+\zeta(7).
\end{align*}
Now due to the algebraic relations between multiple zeta values (the shuffle and stuffle relations), all multiple zeta values of weight $7$ or less can be computed as a weight preserving $\mathbb{Q}$-linear combination of products of single zeta values. This follows from writing out the relations found in theorems 3.1, 3.2, 3.3 in these lecture notes by Wadim Zudilin. (The weight of the product $\zeta(s_1)\dots\zeta(s_k)$ is the sum $s_1+\dots+s_k$.)
An advantage of this method is that it works in high generality. For example, if one has a series of the form
$$\sum_{n=1}^\infty\frac{H_n^{(i_1)}H_n^{(i_2)}\ldots H_n^{(i_k)}}{n^s},$$
with $s, i_1,\dots, i_k$ positive integers and $s>1$, then it can be written as a $\mathbb{Z}$-linear combination of multiple zeta values of weight $w=s+i_1+\dots+i_k$. Therefore, if $w\leq 7$, then the series can be written as a $\mathbb{Q}$-linear combination of products of single zeta values of weight $w$.
