Compute: $\sum_{k=1}^{\infty}\sum_{n=1}^{\infty} \frac{1}{k^2n+2nk+n^2k}$

I try to solve the following sum:

$$\sum_{k=1}^{\infty}\sum_{n=1}^{\infty} \frac{1}{k^2n+2nk+n^2k}$$

I'm very curious about the possible approaching ways that lead us to solve it. I'm not experienced with these sums, and any hint, suggestion is very welcome. Thanks.

I think one way of approaching this sum would be to use the partial fraction $$ \frac{1}{k^2n+2nk+n^2k} = \frac{1}{kn(k+n+2)} = \frac{1}{2}\Big(\frac{1}{k} + \frac{1}{n}\Big)\Big(\frac{1}{k+n} - \frac{1}{n+k+2}\Big)$$ to rewrite you sum in the form $$\sum_{n=1}^{\infty}\sum_{k=1}^{\infty} \frac{1}{k^2n+n^2k+2kn} = \frac{1}{2}\sum_{k=1}^{\infty}\sum_{n=1}^{\infty} \Big( \frac{1}{n(k+n)} - \frac{1}{n(k+n+2)} + \frac{1}{k(k+n)} - \frac{1}{k(k+n+2)}\Big)$$ Since the sum on the right will telescope in one of the summation variables it should be straightforward to find the answer from here (it ends up being $7/4$ I think).

Let's write: $$f(x) = \sum_{k = 1}^{+\infty}\sum_{n = 1}^{+\infty} \frac{x^{n+k+2}}{nk (n+k+2)}$$ then: $$f'(x) = \sum_{k = 1}^{+\infty}\sum_{n = 1}^{+\infty} \frac{x^{n+k+1}}{nk} = - \sum_{k = 1}^{+\infty}\frac{x^{k+1} \ln (1-x)}{k}$$ We want to know the value of $f(1)$, so we integrate: $$f(1) = - \sum_{k = 1}^{+\infty}\frac{1}{k} \int_0^1 x^{k+1} \ln (1-x) \, dx = \sum_{k = 1}^{+\infty} \frac{H_{k+2}}{k(k+2)}\\ = \frac{1}{2} \sum_{k = 1}^{+\infty} \left( \frac{H_{k+2}}{k} - \frac{H_{k+2}}{k+2} \right)=\frac{1}{2} \sum_{k = 1}^{+\infty} \left( \frac{H_{k} + \frac{1}{k+1} + \frac{1}{k+2}}{k} - \frac{H_{k+2}}{k+2} \right) \\ = \frac{1}{2} \left( H_1 + \frac{1}{2}H_2 + \sum_{k = 1}^{+\infty} \left( \frac{3}{2k} - \frac{1}{1+k} - \frac{1}{2(k+2)} \right) \right) = \frac{7}{4}$$ The very last sum telescopes: $$\sum_{k = 1}^{+\infty} \left( \frac{3}{2k} - \frac{1}{1+k} - \frac{1}{2(k+2)} \right) = \sum_{k = 1}^{+\infty} \left( \left( \frac{1}{k} - \frac{1}{1+k} \right) + \left(\frac{1}{2k} - \frac{1}{2(k+2)} \right) \right) = \frac{7}{4}$$ And the integral I've used is obtained by integrating by parts: $$\int_0^1 x^{k+1} \ln(1-x)\,dx = \frac{x^{k+2} \ln(1-x)}{k+2} \Big|_0^1 + \frac{1}{k+2} \int_0^1 \frac{x^{k+2} - 1 + 1}{1-x} \,dx \\ = -\frac{H_{k+2}}{k+2} + \lim_{\epsilon \to 1^-} \left( x^{k+2} \ln(1-x) + \int_0^\epsilon \frac{dx}{1-x} \right) = -\frac{H_{k+2}}{k+2}$$

I admit it is a little bit long solution but the trick of making power series from regular series is very useful :)

Here's another approach. It depends primarily on the properties of telescoping series, partial fraction expansion, and the following identity for the $m$th harmonic number $$\begin{eqnarray*} \sum_{k=1}^\infty \frac{1}{k(k+m)} &=& \frac{1}{m}\sum_{k=1}^\infty \left(\frac{1}{k} - \frac{1}{k+m}\right) \\ &=& \frac{1}{m}\sum_{k=1}^m \frac{1}{k} \\ &=& \frac{H_m}{m}, \end{eqnarray*}$$ where $m=1,2,\ldots$.

Then,
$$\begin{eqnarray*} \sum_{k=1}^{\infty}\sum_{n=1}^{\infty} \frac{1}{k^2n+2nk+n^2k} &=& \sum_{k=1}^{\infty} \frac{1}{k} \sum_{n=1}^{\infty} \frac{1}{n(n+k+2)} \\ &=& \sum_{k=1}^{\infty} \frac{1}{k} \frac{H_{k+2}}{k+2} \\ &=& \frac{1}{2} \sum_{k=1}^{\infty} \left( \frac{H_{k+2}}{k} - \frac{H_{k+2}}{k+2} \right) \\ &=& \frac{1}{2} \sum_{k=1}^{\infty} \left( \frac{H_k +\frac{1}{k+1}+\frac{1}{k+2}}{k} - \frac{H_{k+2}}{k+2} \right) \\ &=& \frac{1}{2} \sum_{k=1}^{\infty} \left( \frac{H_k}{k} - \frac{H_{k+2}}{k+2} \right) + \frac{1}{2} \sum_{k=1}^{\infty} \left(\frac{1}{k(k+1)} + \frac{1}{k(k+2)}\right) \\ &=& \frac{1}{2}\left(H_1 + \frac{H_2}{2}\right) + \frac{1}{2}\left(H_1 + \frac{H_2}{2}\right) \\ &=& \frac{7}{4}. \end{eqnarray*}$$