I was trying to get some interesting result for $\zeta(3)$, exploring the following function:
$$W(a) = \sum_{k=1}^\infty \frac{1}{k^3 + a^3}, \mbox{ with } \lim_{a\rightarrow 0} W(a) = \zeta(3).$$

Let $w_1, w_2, w_3$ be the three roots (one real, two complex) of $(w+1)^3+a^3=0$, with $w_1=-(a+1)$. Also, $a$ is a real number. Using Wolfram Alpha (see [here][1]), I get

$$W(a)=\frac{-1}{3}\cdot\sum_{j=1}^3 W_j(a), \mbox{ with } W_j(a) = \frac{\psi^{(0)}(-w_j)}{(w_j+1)^2}.$$

Here $\psi^{(0)}$ is the digamma function. The result is wrong because $W_1(a) \rightarrow \infty$ as $a\rightarrow 0^+$ while $W_2(a)$ and $W_3(a)$ remain bounded. Indeed using $a=0.0001$, Wolfram yields $W(a)\approx -2334.16$, see [here][2]. Surprising, with $a=0.01$ it yields $W(a)\approx 1.20206$ which is very close to the true result. 

Surprisingly, Wolfram knows (see [here][3]) that 
$$\lim_{a\rightarrow 0} W(a) = -\frac{\psi^{(2)}(1)}{2}.$$

Of course (this is a well known fact), $\zeta(3)=-\psi^{(2)}(1)/2$ and thus Wolfram is correct this time.

**My question:** 

What is going on with this computation (or is it me?), and what is the correct formula for $W(a)$?

**Update**

See the two answers below proving that I was wrong, and that the Mathematica formula I though was incorrect, is indeed right. Kudos Mathematica! You were successful at solving a nice problem involving a few challenging steps, and coming with a somewhat unexpected but neat formula involving derivatives of the digamma function instead of the classic $\zeta(3)$. 

**Final note**

It is possible to use a different, simpler approach that does not involve complex numbers. Consider

$$V(a) =\sum_{k=1}^\infty (-1)^{k+1}\frac{1}{k(k^2-a^2)}.$$

Wolfram is able to compute the limit of $V(a)$ as $a\rightarrow 0$, and returns the correct value $3\zeta(3)/4$, see [here][4]. It is easy to establish that

$$V(a)=\frac{1}{a^2} \Big[\int_0^\infty \frac{\cosh(ax)}{1+e^x} dx -\log 2\Big].$$

To compute $\lim_{a\rightarrow 0} V(a)$, we apply L'Hospital Rule twice to the above expression, the denominator in this case being $a^2$. This yields

$$\lim_{a\rightarrow 0}V(a) = \frac{1}{2}\lim_{a\rightarrow 0}\int_0^\infty \frac{x^2\cosh(ax)}{1+e^x}dx =\frac{1}{2}\int_0^\infty \frac{x^2}{1+e^x}dx=\frac{3\zeta(3)}{4}.$$

Here, we assume $a<1$.


  [1]: https://www.wolframalpha.com/input/?i=sum%201%2F%28k%5E3%20%2Ba%5E3%29%20for%20k%3D1..infinity
  [2]: https://www.wolframalpha.com/input/?i=sum%201%2F%28k%5E3%20%2B%280.0001%29%5E3%29%20for%20k%3D1..infinity
  [3]: https://www.wolframalpha.com/input/?i=lim%28sum%201%2F%28k%5E3%20%2B%20a%5E3%29%20for%20k%3D1..infinity%29%2C%20a%3D0
  [4]: https://www.wolframalpha.com/input/?i=lim%20%28sum%20%28-1%29%5E%28k%2B1%29%2F%28k%5E3-a%5E2%20k%29%2C%20k%3D1..infinity%29%2C%20a%20%3D%200