$\sum_{-\infty}^{+\infty}\frac{\exp(-n^2)}{1-4n^2}$ in closed form?
EDIT 29.06.18 13:15
Correction
The previously provided closed expression $s_{c}$ was shown to be incorrect by some commenters.
As Mariusz pointed out correctly there are jumps in the function $p(z)$ (see (d8)) at $z=k \;\pi$
Correcting for $m$ of those jumps the closed expression is given by the (extremely fast converging) series
$$s_{c}(m) = s_{c}+\frac{\pi }{2 \sqrt[4]{e}} \sum _{k=1}^m (-1)^k \left(\text{erfi}\left(\frac{1}{2}+i \pi k\right)+\text{erfi}\left(\frac{1}{2}-i \pi k\right)\right)$$
Numerical examples
Letting $sN(100) = s$ up to 100 valid digits we obtain for the differences $d(m) = s_{c}(m) - sN(100)$ for $m=0..3$ the following
$$ \left\{-\frac{3.996}{10^6},-\frac{1.54}{10^{19}},-\frac{2.5929}{10^{41}},-\frac{1.456}{10^{71}}\right\}$$
Summarizing: instead of a true closed expression for $s$ we have effectively transformed one sum into another sum, which, however converges very fast.
The original question is now to be repeated for the latter sum.
Remark
Independently of this correction, the discovered strange numerical behaviour Mathematica deserves further study.
Note added 1. July: Some steps have already been taken here https://mathematica.stackexchange.com/questions/176240/bug-in-analytical-expression-of-integral-containing-abs-function/176342#176342
EDIT 28.06.18 14:00
Doubts have been raised in comments that my result might be wrong. The arguments given seem to rely on the numerical evaluation in Mathematica. Therefore I have extended the derivation to show more steps so that possible flaws can be detected.
Post as of 27.06.18
I have found the closed expression for the sum
$$s=\sum _{n=-\infty }^{\infty } \frac{\exp \left(-n^2\right)}{1-4 n^2}$$
It is given by
$$s_{c}=\frac{\pi\; \text{erfi}\left(\frac{1}{2}\right)}{2 \sqrt[4]{e}} $$
Derivation see below.
Numerically, the first 60 digits of $s_{c}$ according to Mathematica 10.1 are
$$N(s_{c}) = 0.75229|3902402569849043685417920199342618157039554017947019766$$
while, as has been pointed out in two comments, $s_{c}$ differs from $s$ numerically already in the fifth decimal digit:
$$N(s) = 0.75229|7898472243144830594123416216568518483359108518774883675$$
Something is wrong here.
Original post
As a partly solution we give here an integral representation of the sum in question
$$s=\sum _{n=-\infty }^{\infty } \frac{\exp \left(-n^2\right)}{1-4 n^2}$$
Splitting the sum into the term with $n=0$ and observing the symmetry of the summands $s$ can be witten as
$$s = 1-2 p\tag{1}$$
where
$$p = \sum _{n=1}^{\infty } \frac{\exp \left(-n^2\right)}{4 n^2-1}\tag{2}$$
Writing the denomintor for $n \ge 1$ as an integral
$$\frac{1}{4 n^2-1} = \int_0^{\infty } \exp \left(-t \left(4 n^2-1\right) \right) \, dt \tag{3}$$
leads under the integral to the sum
$$ \sum _{n=1}^{\infty } \exp \left(-\left(4 n^2-1\right) t-n^2\right) $$
Which can be evaluated in terms of the Jacobi Theta function:
$$\frac{1}{2} e^t \left(\vartheta _3\left(0,e^{-4 t-1}\right)-1\right)\tag{4}$$
Using (3) we find the integral representation
$$p=\int_0^{\infty } \frac{1}{2} e^t \left(\vartheta _3\left(0,e^{-4 t-1}\right)-1\right) \, dt\tag{5}$$
Derivation of the closed expression
The ingredients are easy to verify:
Partial fraction decomposition gives instead of (3):
$$\frac{1}{4 n^2-1} = \int_0^{\infty } \sinh (t) \exp (-2 n t) \, dt\tag{d1}$$
Under the Fourier transform the exponent becomes linear in $n$
$$e^{-n^2} = \frac{1}{\sqrt{\pi }} \int_{-\infty }^{\infty } \exp \left(-x^2\right) \exp (2 i n x) \, dx\tag{d2}$$
Under the two integrals the sum for $p$ is linear in the Exponent, i.e. it is a geometric sum, and hence can easily be done:
$$\frac{1}{\sqrt{\pi }}\sum _{n=1}^{\infty } \exp \left(-x^2\right) \sinh (t) \exp (-2 n t) \exp (2 i n x) = \frac{e^{-x^2+2 i x} \sinh (t)}{\sqrt{\pi } \left(e^{2 t}-e^{2 i x}\right)}\tag{d3}$$
Now the t-integral is solved by Matematica assuming that $x$ is real to give
$$\int_0^{\infty } \frac{e^{-x^2+2 i x} \sinh (t)}{\sqrt{\pi } \left(e^{2 t}-e^{2 i x}\right)} \, dt = \frac{e^{-x^2} \left(1+2 i \sin (x) \tanh ^{-1}\left(e^{i x}\right)\right)}{2 \sqrt{\pi }}\tag{d4}$$
Last but not least, the x-integral between $-\infty$ and $+\infty$ can be written, using the symmetry of the integrand as
$$p=\int_0^{\infty } \frac{e^{-x^2} \left(1-i \sin (x) \left(\tanh ^{-1}\left(e^{-i x}\right)-\tanh ^{-1}\left(e^{i x}\right)\right)\right)}{\sqrt{\pi }} \, dx\tag{d5}$$
Observing that
$$i \left(\tanh ^{-1}\left(e^{-i x}\right)-\tanh ^{-1}\left(e^{i x}\right)\right)=\frac{\pi }{2} \; \text{sgn}(\sin (x))\tag{d6}$$
the x-integral becomes
$$p=\frac{1}{2 \sqrt{\pi }}\int_{-\infty }^{\infty } e^{-x^2} \left(1-\frac{1}{2} \pi \left| \sin (x)\right| \right) \, dx\tag{d7}$$
where we have returned to the symmetric form, correcting this by the factor $\frac{1}{2}$ in front
Mathematica refused to do this integral directly but it was successful with the finite integral, with some $z\gt 0$
$$p(z)=\frac{1}{2 \sqrt{\pi }}\int_{-z }^{z} e^{-x^2} \left(1-\frac{1}{2} \pi \left| \sin (x)\right| \right) \, dx$$
$$= \frac{1}{8} \left(4 \text{erf}(z)+\frac{\pi \left(-2 \text{erfi}\left(\frac{1}{2}\right)+\left(\text{erfi}\left(\frac{1}{2}+i z\right)+\text{erfi}\left(\frac{1}{2}-i z\right)\right) \csc (z) \left| \sin (z)\right| \right)}{\sqrt[4]{e}}\right)\tag{d8}$$
Now the Limit $z\to\infty$ has to be taken. Again Mathematica refused but it did this asymtotic series expansion about $z=\infty$
$$\frac{i \sqrt{\pi } e^{-z^2-i z} \csc (z) \left| \sin (z)\right| }{8 z}-\frac{i \sqrt{\pi } e^{-z^2+i z} \csc (z) \left| \sin (z)\right| }{8 z}-\frac{\pi \text{erfi}\left(\frac{1}{2}\right)}{4 \sqrt[4]{e}}-\frac{e^{-z^2}}{2 \sqrt{\pi } z}+\frac{1}{2}\tag{d9}$$
Letting now $z\to\infty$ gives finally
$$p = \frac{1}{2}-\frac{\pi \text{erfi}\left(\frac{1}{2}\right)}{4 \sqrt[4]{e}}\tag{d10}$$
and the intial sum becomes
$$s = 1 - 2 p = \frac{\pi \text{erfi}\left(\frac{1}{2}\right)}{2 \sqrt[4]{e}}$$
I found general formula. It's not closed form solution only approximation for the sum:
$$\sum _{j=-\infty }^{\infty } \frac{\exp \left(-j^2\right)}{1-4 j^2}\approx \frac{\pi \sum _{k=1}^n \left(\frac{\text{erfi}\left(\frac{1}{2}\right)}{n}+(-1)^k \text{erfi}\left(\frac{1}{2}-k i \pi \right)+(-1)^k \text{erfi}\left(\frac{1}{2}+k i \pi \right)\right)}{2 \sqrt[4]{e}}$$
for $n=1$ it's only correct for 18 digits.
for $n=2$ it's only correct for 40 digits.
for $n=3$ it's only correct for 70 digits.
for $n=4$ it's only correct for 109 digits.
if $n>4$ then it's better approximation and more correct digits.
Derivation of the formula
Borrowing integral $(d7)$ form user Dr. Wolfgang Hintze and integrating with Mathematica help:
$$\int \frac{\exp \left(-x^2\right) \left(1-\frac{1}{2} \pi \left| \sin (x)\right| \right)}{2 \sqrt{\pi }} \, dx=\\\frac{4 \sqrt[4]{e} \text{erf}(x)-\pi \left(\text{erfi}\left(\frac{1}{2}-i x\right)+\text{erfi}\left(\frac{1}{2}+i x\right)\right)+2 \pi \left(\text{erfi}\left(\frac{1}{2}-i x\right)+\text{erfi}\left(\frac{1}{2}+i x\right)\right) \theta (\sin (x))}{16 \sqrt[4]{e}}+C$$
taking limit in jumps points:
$$\left(\underset{x\to \pi ^-}{\text{lim}}\text{int}-\underset{x\to 0^+}{\text{lim}}\text{int}\right)+\left(\underset{x\to \infty }{\text{lim}}\text{int}-\underset{x\to (2 \pi )^+}{\text{lim}}\text{int}\right)+\left(\underset{x\to (2 \pi )^-}{\text{lim}}\text{int}-\underset{x\to \pi ^+}{\text{lim}}\text{int}\right)+\left(\underset{x\to (-2 \pi )^-}{\text{lim}}\text{int}-\underset{x\to -\infty }{\text{lim}}\text{int}\right)+\left(\underset{x\to 0^-}{\text{lim}}\text{int}-\underset{x\to (-\pi )^+}{\text{lim}}\text{int}\right)+\left(\underset{x\to (-\pi )^-}{\text{lim}}\text{int}-\underset{x\to (-2 \pi )^+}{\text{lim}}\text{int}\right)=\\\frac{1}{2}-\frac{\pi \left(\text{erfi}\left(\frac{1}{2}\right)-\text{erfi}\left(\frac{1}{2}-i \pi \right)-\text{erfi}\left(\frac{1}{2}+i \pi \right)+\text{erfi}\left(\frac{1}{2}-2 i \pi \right)+\text{erfi}\left(\frac{1}{2}+2 i \pi \right)\right)}{4 \sqrt[4]{e}}$$
based on this, it was possible to deduce the formula.