$$S=\sum_{r=1}^{\infty}\sum_{s=1}^{\infty} \frac{1}{rs(r+s)}=\int_{0}^1 \sum_{r=1}^{\infty}\sum_{s=1}^{\infty}\frac{x^{r+s-1}}{rs} dx$$ $$=\int_{0}^{1}\frac{\log^2(1-x)}{x} dx= \int_{0}^{1} \frac{\log^2 t}{1-t}dt= \int_{0}^{\infty} \frac{u^2}{e^u-1} du=2\zeta(3)=2\sum_{k=1}^{\infty}\frac{1}{k^3}.$$
The question is: what are other ways of proving the above summation without using integration?
We group terms in the double sum by setting $n=r+s$. For each $n\geq2$, the pairs of positive integers $(r,s)$ satisfying $r+s=n$ are exactly those with $r=1$ to $n−1$, and $s=n−r$. With this we rewrite the double sum as a single sum over $n$ with an inner sum over $r$.
$$ \sum_{r=1}^\infty \sum_{s=1}^\infty \frac{1}{rs(r+s)} = \sum_{n=2}^\infty \sum_{r=1}^{n-1} \frac{1}{r(n - r)n} $$
We can then manipulate the expression
$$ \frac{1}{r(n - r)} = \frac{1}{n} \left( \frac{1}{r} + \frac{1}{n - r} \right) $$
So we get
$$ \sum_{r=1}^{n-1} \frac{1}{r(n - r)} = \frac{1}{n} \sum_{r=1}^{n-1} \left( \frac{1}{r} + \frac{1}{n - r} \right) = \frac{2}{n} \sum_{r=1}^{n-1} \frac{1}{r} = \frac{2H_{n-1}}{n} $$
Therefore we can write that
$$ \sum_{r=1}^\infty \sum_{s=1}^\infty \frac{1}{rs(r+s)} = \sum_{n=2}^\infty \frac{1}{n} \cdot \frac{2H_{n-1}}{n} = 2 \sum_{n=2}^\infty \frac{H_{n-1}}{n^2} $$
We now shift the index by letting $m = n - 1$, so:
$$ 2 \sum_{n=2}^\infty \frac{H_{n-1}}{n^2} = 2 \sum_{m=1}^\infty \frac{H_m}{(m+1)^2} $$
It is known(1.12, page 3) that
$$ \sum_{m=1}^\infty \frac{H_m}{(m+1)^2} = \zeta(3) $$
Hence, we can conclude that
$$ \sum_{r=1}^{\infty} \sum_{s=1}^{\infty} \frac{1}{rs(r+s)} = 2\zeta(3) $$
