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

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

I came across these two sums while working on an tough one of wight 9 and I managed to find these two results but I don't think my solution is a good approach as it's pretty long and involves tedious calculations, so I am seeking different methods if possible. I am much into new ideas. All approaches are appreciated though. 
By the way, do we have a generalization for $\displaystyle\sum_{n=1}^\infty \frac{H_n^{(2)}}{n^a}$, for odd $a$? Note that  there is no closed form for even $a>4$.
Thanks in advance.
Note: You can find these two results on Wolfram Alpha here and here respectively but I modified it a little bit as I like it expressed in terms of $\zeta(a)$ instead of $\pi^a$.
 A: Your method can be carried out in full generality as long as the power $a$ is odd. Starting with
$$\sum_{n=1}^\infty \frac{H_n^{(2)}}{n^a} = \zeta(2)\zeta(a)-\sum_{n,m=1}^\infty \frac{1}{n^a(n+m)^2} \, ,$$ we can generally decompose the fraction in the second term into partial fractions $$\frac{1}{n^a(n+m)^2}=\sum_{k=1}^a \frac{A_k}{n^k}+\frac{B_1}{n+m} + \frac{B_2}{(n+m)^2} \, . \tag{1}$$
By the method of residues these coefficients are given by $$A_k=\frac{(-1)^{a-k}(a-k+1)}{m^{a-k+2}} \\ B_1=\frac{a\,(-1)^a}{m^{a+1}} \\ B_2=\frac{(-1)^a}{m^a}$$
and hence
$$\sum_{n=1}^\infty \frac{H_n^{(2)}}{n^a} = \zeta(2)\zeta(a)-\sum_{n,m=1}^\infty \frac{a \, (-1)^{a-1}}{m^{a+1}} \left(\frac{1}{n} - \frac{1}{m+n} \right) \\
-\sum_{n,m=1}^\infty \sum_{k=2}^a \frac{(-1)^{a-k}(a-k+1)}{n^k m^{a-k+2}} - \sum_{n,m=1}^\infty \frac{(-1)^a}{m^a (m+n)^2} \\
= \zeta(2)\zeta(a) + a \, (-1)^a \sum_{m=1}^\infty \frac{H_m}{m^{a+1}} \\
- \sum_{k=2}^a (-1)^{a-k} (a-k+1) \zeta(k)\zeta(a-k+2) - (-1)^a \left(\zeta(2)\zeta(a) - \sum_{m=1}^\infty \frac{H_m^{(2)}}{m^a} \right) \, .$$ For even $a$ you will get an identity similar to that of Euler. For odd $a$ you can solve for the LHS $$\sum_{n=1}^\infty \frac{H_n^{(2)}}{n^a} = \zeta(2)\zeta(a) - \frac{a}{2} \sum_{m=1}^\infty \frac{H_m}{m^{a+1}} + \frac{1}{2} \sum_{k=2}^a (-1)^{k} (a-k+1) \zeta(k)\zeta(a-k+2) \, .$$
You can simplify it further by substituting Euler's identity for the middle term.

The coefficients in (1) are obtained as follows. For $A_k$ multiply (1) by $n^a$ and derive both sides with respect to $n$ exactly $a-k$ times. Finally set $n=0$. This way the LHS will yield $\frac{(-1)^{a-k} \, (a-k+1)!}{m^{a-k+2}}$. The $B$-terms on the RHS will vanish, because $n^a$ is derived at most $a-k<a$ times, which leaves at least one power in $n$ which is set to zero. For the $A$-terms the monomials $n^{a-k}$ are successively derived and those with power $<a-k$ vanish immediately, while those with power $>a-k$ will vanish upon setting $n=0$. The term with power $a-k$ gives $(a-k)! A_k$.
The $B$ coefficients are obtained in the same way by multiplying (1) with $(n+m)^2$, deriving zero or one times and setting $n=-m$.
