<s>[INST] A series involving the harmonic numbers : $\sum_{n=1}^{\infty}\frac{H_n}{n^3}$ Let $H_{n}$ be the nth harmonic number defined by $ H_{n} := \sum_{k=1}^{n} \frac{1}{k}$.
How would you prove that 
$$\sum_{n=1}^{\infty}\frac{H_n}{n^3}=\frac{\pi^4}{72}?$$

Simply replacing $H_{n}$ with $\sum_{k=1}^{n} \frac{1}{k}$ does not seem like a good starting point. Perhaps another representation of the nth harmonic number would be more useful.
 [/INST] I will try to reduce the sum to an integral:
$$
 \sum_{n=1}^\infty \frac{H_n}{n^3} = \sum_{n=1}^\infty H_n \frac{1}{\Gamma(3)} \int_0^\infty x^2 \mathrm{e}^{-n x} \mathrm{d} x = \frac{1}{2} \int_0^\infty x^2 \sum_{n=1}^\infty H_n \mathrm{e}^{-n x} \mathrm{d} x \tag{1}
$$
We now make use of the following generating function:
$$
   \sum_{n=1}^\infty H_n z^n = \sum_{n=1}^\infty H_n \Delta_n \left(\frac{z^n}{z-1} \right)
$$
where $\Delta_n f_n = f_{n+1}-f_n$. We can now use summation by parts:
$$
     \sum_{n=1}^m a_n \Delta_n b_n = b_{m+1} a_m - b_1 a_1 - \sum_{n=1}^{m-1} b_{n+1} \Delta_n a_n
$$
with $b_n = \frac{z^n}{z-1}$ and $a_n = H_n$, and using $\Delta_n H_n = \frac{1}{n+1}$, we get
$$
  \sum_{n=1}^\infty H_n z^n =  \sum_{n=1}^\infty H_n \Delta_n \left(\frac{z^n}{z-1} \right) = -1 - \sum_{n=1}^\infty \frac{z^{n+1}}{z-1} \frac{1}{n+1} = \frac{\log(1-z)}{z-1}
\tag{2} $$
Now, using $(2)$ in $(1)$:
$$
  \sum_{n=1}^\infty \frac{H_n}{n^3} = -\frac{1}{2} \int_0^\infty x^2 \frac{\log\left(1-\mathrm{e}^{-x}\right)}{1-\mathrm{e}^{-x}} \mathrm{d}x \stackrel{t=\exp(-x)}{=} -\frac{1}{2} \int_0^1 \frac{\log(1-t)}{1-t} \frac{\log^2(t)}{t} \mathrm{d}t \tag{3}
$$ 
The latter integral can be evaluated using derivatives of the Euler beta function:
$$
   \int_0^1 \frac{\log(1-t)}{1-t} \frac{\log^2(t)}{t} \mathrm{d}t = \lim_{\alpha \downarrow 0} \lim_{\beta \downarrow 0} \frac{\mathrm{d}}{\mathrm{d} \alpha} \frac{\mathrm{d}^2}{\mathrm{d} \beta^2} \int_0^1 \left(1-t\right)^{\alpha-1} t^{\beta-1} \mathrm{d} t = \lim_{\alpha \downarrow 0} \lim_{\beta \downarrow 0} \frac{\mathrm{d}}{\mathrm{d} \alpha} \frac{\mathrm{d}^2}{\mathrm{d} \beta^2}  \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha + \beta)}
$$
Using 
$$
  \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha + \beta)} = \left(\frac{1}{\alpha} + \frac{1}{\beta} \right) \frac{\Gamma(\alpha+1) \Gamma(\beta+1)}{\Gamma(\alpha + \beta+1)} = \left(\frac{1}{\alpha} + \frac{1}{\beta} \right)  \left( 1 - \frac{\pi^2}{6} \alpha \beta + \left(\alpha \beta^2 + \beta \alpha^2\right) \zeta(3) - \frac{\pi^4}{360} \left(4 \alpha \beta^3 + \alpha^2 \beta^2 + 4 \alpha^3 \beta\right) + \cdots \right)
$$
Differentiating we get the result:
$$
    \lim_{\alpha \downarrow 0} \lim_{\beta \downarrow 0} \frac{\mathrm{d}}{\mathrm{d} \alpha} \frac{\mathrm{d}^2}{\mathrm{d} \beta^2}  \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha + \beta)} = -\frac{\pi^2}{36}
$$
yielding with eq. $(3)$:
$$
   \sum_{n=1}^\infty \frac{H_n}{n^3} = \frac{\pi^4}{72}
$$
</s>