By using the integral representation of harmonic numbers and by using elementary integration we can fairly easily find the generating function of squares of harmonic numbers: \begin{equation} {\mathcal S}(x) := \sum\limits_{n=1}^\infty H_n^2 x^n = \frac{Li_2(1)- Li_2(1-x)}{1-x} + \frac{-\log(x) \log(1-x) + [\log(1-x)]^2}{1-x} \end{equation} Now let us take $q \ge 2$. By consecutively dividing both sides by $x$ and then integrating $q$ times we easily obtain the following identity: \begin{eqnarray} &&(-1)^{q-1} (q-1)!\sum\limits_{n=1}^\infty \frac{H_n^2}{n^q} = \\ &&\int\limits_0^1 \left(\frac{Li_2(1)- Li_2(1-x)}{x}\right) \frac{[\log(x)]^{q-1}}{1-x} dx +\int\limits_0^1 \frac{\left(-[\log(x)]^q \log(1-x) + [\log(x)]^{q-1} [\log(1-x)]^2\right)}{x(1-x)} dx \end{eqnarray} The second integral on the right hand side is easily expressed through integrals given in Compute an integral containing a product of powers of logarithms. and is a function of polygamma function values at unity. Now the question would be how do we compute the first integral? By using the definition of the dilogarithm and by integrating by parts it seems plausible to do that integral for $q=2$ or maybe $q=3$ and therefore find the Euler sum in question. Can we compute that integral for generic values of $q$ ?
By using the reflection formula for the dilogarithmic function we simplify the right hand side as follows: \begin{eqnarray} rhs &=& \int\limits_0^1 \frac{Li_2(x)}{x} \cdot \frac{[\log(x)]^{q-1}}{(1-x)} dx &+& \int\limits_0^1 \frac{[\log(x)]^{q-1} [\log(1-x)]^2}{x(1-x)} dx \\ &=& \int\limits_0^1 Li_2(x) \cdot [\log(x)]^{q-1}\cdot \left( \frac{1}{x} + \frac{1}{1-x}\right) dx &+& \left({\mathcal I}^{(2,q-1)} + {\mathcal I}^{(q-1,2)}\right) \\ &=& (-1)^{q-1}(q-1)! \left(Li_{q+2}(1)-Li_{q+2}(1)+\sum\limits_{n=1}^\infty \frac{H_n^{(2)}}{n^q}\right) &+&\left({\mathcal I}^{(2,q-1)} + {\mathcal I}^{(q-1,2)}\right) \\ \end{eqnarray} In the second line we introduced: \begin{equation} {\mathcal I}^{(q,p)} := \int\limits_{0}^1 \frac{[\log(1-x)]^q [\log(x)]^p}{x} dx \end{equation} and in the last line we expanded the terms $Li_2(x)/x$ and $Li_2(x)/(1-x)$ in a series in $x$ and integrated term by term. Therefore we get: \begin{equation} \sum\limits_{n=1}^\infty \frac{H_n^2}{n^q} = \sum\limits_{n=1}^\infty \frac{H_n^{(2)}}{n^q} + \frac{(-1)^{q-1}}{(q-1)!} \cdot \left({\mathcal I}^{(2,q-1)} + {\mathcal I}^{(q-1,2)}\right) \end{equation} The quantities ${\mathcal I}^{(q,p)}$ have been calculated in Compute an integral containing a product of powers of logarithms. . We have: \begin{eqnarray} &&\sum\limits_{m=1}^\infty \frac{[H_m]^2 - H_m^{(2)}}{m^q} = \\ &&\frac{(-1)^{q-1}}{(q-1)!} \left(\right.\\ && -(\frac{1}{q}+\frac{1}{3}) \Psi^{(q+1)}(1) + \sum\limits_{j=1}^{q-1} (\frac{1}{q} \binom{q}{j} + \frac{1_{q\ge 3}}{2}(\binom{q}{j}-1_{j\notin[2,q-2]})) \Psi^{(j)}(1) \Psi^{(q-j)}(1)-\\ &&\frac{1}{3} \sum\limits_{1\le j < j_1 \le q-2} \frac{(q-1)!}{j!(j_1-j)!(q-1-j_1)!} \Psi^{(j)}(1) \Psi^{(j_1-j)}(1) \Psi^{(q-1-j_1)}(1)\\ && \left.\right) \end{eqnarray} for $q\ge 2$. Particular cases of the sum above are given in http://algo.inria.fr/flajolet/Publications/FlSa98.pdf .
