I am trying to use the Euler-Maclaurin formula to obtain the following difference:
\begin{align} \Delta &= \sum_{n=1}^\infty f(n) - \int_0^\infty dx f(x) \\ &= \frac{ f(\infty) - f(0) }{2} + \sum_{k=1}^{\lfloor{p/2}\rfloor} \frac{ B_{2k}}{(2k)!} \left( f^{2k-1}(\infty) - f^{2k-1}(0) \right) + R_p \end{align} with \begin{align} f(x) = \frac{1}{(x^2 + a^2)^2}. \end{align}
$f$ and all its derivatives of this function are zero at $x\to\infty$ and, moreover, since it is an even function, it's odd derivatives vanish at $x=0$. On the other hand, from the bound on $R_p$ given in Wikipedia, it seems $R_{p \to \infty} = 0$. Therefore, I get \begin{align} \Delta = -\frac{f(0)}{2} = -\frac{1}{2 a^4} \end{align}
However, Mathematica gives me exact expression for the sum, and the integral can also be carried out \begin{align} \sum_{n=1}^\infty f(n) &= \frac{\pi a \left(\coth (\pi a)+\pi a \, \text{csch}^2(\pi a)\right)-2}{4 a^4}, \qquad\qquad \int_0^\infty dx \, f(x) = \frac{\pi}{4 a^3}, \end{align} for which the difference seems not to agree with $\Delta$ calculated using Euler-Maclaurin. What am I missing or doing wrong here?