Correct method of integration involving two exponential terms I have an Integrand involving two exponential terms:
$$
\int_{0}^{\infty} \frac{\exp(x^2)}{(1+\exp(x^2))^2} dx
$$
I would like to know what is the best way to integrate such a function without blowing it up? 
What if $x^2$ is replaced by two variables $(x^2 + y^2)$ and we have a double integral? Will the method of integration remain the same?
I use python and matlab for calculations.
Thanks
Notation fixed.

$$
\begin{align}
\int_0^\infty\frac{e^{x^2}}{\left(1+e^{x^2}\right)^2}\,\mathrm{d}x
&=\int_0^\infty\frac{e^{-x^2}}{\left(1+e^{-x^2}\right)^2}\,\mathrm{d}x\\
&=\frac12\int_0^\infty\frac{e^{-x}}{(1+e^{-x})^2}\frac{\mathrm{d}x}{\sqrt{x}}\tag{1}
\end{align}
$$
Consider
$$
\begin{align}
&\int_0^\infty\frac{e^{-x}}{(1+e^{-x})^2}x^{\alpha-1}\mathrm{d}x\tag{2}\\
&=\int_0^\infty\left(e^{-x}-2e^{-2x}+3e^{-3x}-\dots\right)x^{\alpha-1}\,\mathrm{d}x\tag{3}\\
&=\Gamma(\alpha)\left(1-2^{1-\alpha}+3^{1-\alpha}-\dots\right)\tag{4}\\[6pt]
&=\Gamma(\alpha)\zeta(\alpha-1)\left(1-2^{2-\alpha}\right)\tag{5}
\end{align}
$$
The integral in $(2)$ is analytic for $\mathrm{Re}(\alpha)\gt0$. For $\mathrm{Re}(\alpha)\gt1$ the sum in $(4)$ converges and takes the value in $(5)$. Since the functions in $(2)$ and $(5)$ are analytic and equal for $\mathrm{Re}(\alpha)\gt1$, they must be equal for $\mathrm{Re}(\alpha)\gt0$.
Therefore, with $\alpha=\frac12$, we get
$$
\begin{align}
\int_0^\infty\frac{e^{x^2}}{\left(1+e^{x^2}\right)^2}\,\mathrm{d}x
&=\frac{\sqrt\pi}2\zeta\left(-\tfrac12\right)\left(1-\sqrt8\right)\\
&\doteq0.336859119428876991346\tag{6}
\end{align}
$$

A note on computing $\boldsymbol{\zeta\!\left(-\frac12\right)}$
We can't use the standard series
$$
\zeta(s)=\sum_{n=1}^\infty n^{-s}\tag{7}
$$
to compute $\zeta\!\left(-\tfrac12\right)$ because the series in $(7)$ only converges for $s\gt1$.
In this answer, analytic continuation is used to show that 
$$
\zeta\!\left(-\tfrac12\right)=\lim_{n\to\infty}\left(\sum_{k=1}^n\sqrt{k}\,-\tfrac23n^{3/2}-\tfrac12n^{1/2}\right)\tag{8}
$$
The convergence of $(8)$ is very slow; the error is approximately $\frac1{24\sqrt{n}}$. However, by using $8$ terms from the Euler-Maclaurin Sum Formula, the error is reduced to $\frac{52003}{100663296}n^{-25/2}$.
Thus, using $n=1000$, we get
$$
\zeta\!\left(-\tfrac12\right)=-0.2078862249773545660173067253970493022262\dots\tag{9}
$$
