Euler-Maclaurin summation for $e^{-x^2}$ I want to approximate the sum
$$\sum_{k=0}^\infty e^{-k^2}$$
using the Euler-Maclaurin formula
$$\sum_{k=0}^\infty f(k) = \int_0^\infty f(x) \, dx + \frac{1}{2}(f(0) + f(\infty)) + \frac{1}{12}(f'(\infty) - f'(0)) - \frac{1}{720}(f'''(\infty) - f'''(0)) + \ldots$$
where $f(\infty)$ means $\lim_{x \to \infty} f(x)$. But putting $f(x) = e^{-x^2}$, it's not too hard to see that all derivatives of odd order are an odd polynomial times $f$, implying they vanish at $x = 0$. But this means all terms of the Euler-Maclaurin formula vanish except for the first two! But I know this is wrong because I evaluated the sum numerically in Mathematica.
What's wrong with my calculations?
 A: Nothing is wrong with the details of the computation, but you are treating the Euler-Maclaurian formula as an equality, when it isn't. 
Before continuing, I will switch the sum a little bit: You are working with
$$S:=\sum_{k=0}^{\infty} e^{-k^2}.$$
The formulas are a little cleaner in terms of 
$$T:=\sum_{k=-\infty}^{\infty} e^{-k^2}.$$
These are related by $T=2S-1$, so it is easy to switch between them.
The Euler-Macluarin formula with remainder term, combined with your correct computation that $f^{(k)}(x)$ goes to $0$ as $x \to \pm \infty$, tells us that we have
$$\sum_{k=-\infty}^{\infty} e^{-k^2} = \int_{x=-\infty}^{\infty} e^{-x^2} dx + \int_{-\infty}^{\infty} \frac{B_N(x-\lfloor x \rfloor)}{N!}  \frac{d^N e^{-x^2}}{(dx)^N}  dx$$
for any positive integer $N$. Here $B_N$ is the $N$-th Bernouli polynomial and $\lfloor x \rfloor$ is $x$ rounded down to an integer. For some functions, the remainder $\int \frac{B_N(x-\lfloor x \rfloor)}{N!} \frac{d^N f}{(dx)^N} dx$ goes to $0$ as $N \to \infty$, but it often doesn't, and it doesn't in your case. 
To give some other examples, if you compute the asymptotics of $\sum_{k \leq M} \frac{1}{k}$ using Euler-Maclaurin, the Euler-Mascheroni constant occurs as the limit of the remainder term. If you derive Stirling's formula by Euler-Maclaurin summation of $\sum \log k$, then the $\log (2 \pi)$ occurs as the limit of the remainder term.
There is a really good discussion of this in Chapter 9 of Concrete Mathematics, by Graham, Knuth and Patashnik. You might particularly enjoy the fourth example in Chapter 9.6, where they use Euler Maclaurin summation to show that, for any $N$, we have
$$\sum_{k=-\infty}^{\infty} e^{-k^2/t} = \sqrt{\pi t} + O(t^{-N}) \ \mbox{as}\ t \to \infty.$$
A: Many examples of Euler-Maclaurin  summation are actually harmonic sums
and can be treated by Mellin transform methods.
In the present case put
$$S(x) = \sum_{k\ge 1} e^{-x^2 k^2}$$
with so that we are interested in $S(1/\sqrt{t})$
as $t\rightarrow\infty.$
This sum can be evaluated by inverting its Mellin transform.
Recall the harmonic sum identity
$$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x);s\right) =
\left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s)$$
where $g^*(s)$ is the Mellin transform of $g(x).$
In the present case we have
$$\lambda_k = 1, \quad \mu_k = k 
\quad \text{and} \quad
g(x) = e^{-x^2}.$$
We need the Mellin transform $g^*(s)$ of $g(x)$ which is
$$\int_0^\infty e^{-x^2} x^{s-1} dx.$$
Use the substitution $x^2 = u$ so that $2x \; dx = du$ to get
$$\int_0^\infty e^{-u} u^{1/2s-1/2} \frac{1/2 \; du}{\sqrt{u}}
= \frac{1}{2}
\int_0^\infty e^{-u} u^{1/2s-1} du = 
\frac{1}{2} \Gamma(s/2).$$

The  fundamental strip  of this Mellin transform is
$\langle 0,\infty\rangle.$

It  follows that  the Mellin  transform $Q(s)$  of the  harmonic sum
$S(x)$ is given by
$$Q(s) = \frac{1}{2} \Gamma(s/2) \zeta(s)
\quad\text{because}\quad
\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = 
\sum_{k\ge 1} \frac{1}{k^s}
= \zeta(s)$$
for $\Re(s) > 1.$
The Mellin inversion integral for this transform is $$\frac{1}{2\pi i}
\int_{3/2-i\infty}^{3/2+i\infty}  Q(s)/x^s ds$$  which we  evaluate by
shifting it  to the left for  an expansion about zero  (recall that as
$t\rightarrow\infty$ we have $1/\sqrt{t}\rightarrow 0.$)

Observe that the poles are at $s=1$ from the zeta function term and at
the non-positive  even integers from the gamma  function term. However
all of the  latter except the one at zero are  canceled by the trivial
zeros of the  zeta function term, leaving only the  poles at $s=0$ and
$s=1.$

For these two we have
$$\mathrm{Res}\left(Q(s)/x^s; s=1\right) = 
\frac{1}{2} \Gamma(1/2)\frac{1}{x} = \frac{\sqrt{\pi}}{2x}$$
and
$$\mathrm{Res}\left(Q(s)/x^s; s=0\right) =
\frac{1}{2} \times 2 \times  -\frac{1}{2} = -\frac{1}{2}.$$
It follows that as $t\rightarrow\infty$ we have
$$S(1/\sqrt{t}) \sim \frac{1}{2} \sqrt{\pi t} - \frac{1}{2}$$
and in particular
$$2S(1/\sqrt{t})+1 = \sum_{k=-\infty}^\infty e^{-k^2/t}
\sim \sqrt{\pi t}.$$

As for  the error  term if we  have shifted  the integral to  the line
$\Re(s) = -q/2$  with $q>1$ and $q$  odd we have for the  norm of the
zeta function term on the line $-q/2+iv$ the bound $|v|^{1/2+q/2}$ and the
gamma function  term decays  exponentially in $v$ along vertical lines 
and in $q$ at the values  at $-q/2$, so   the   error   term    decays exponentially   also.
