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.
 A: You integral can be written as 
$$I=\frac{1}{4}\int_{0}^{+\infty}\frac{dx}{\cosh^2\left(\frac{x^2}{2}\right)}\,dx=\frac{1}{4\sqrt{2}}\int_{0}^{+\infty}\frac{dz}{\sqrt{z}\cosh^2(z)}=\frac{1}{8\sqrt{2}}\int_{0}^{+\infty}\frac{\tanh(z)}{z^{3/2}}\,dz $$
and due to the identity
$$ \tanh(z) = \sum_{n\geq 0}\frac{8z}{(2n+1)^2 \pi^2+ 4z^2} $$
that follows from considering the logarithmic derivative of the Weierstrass product for the $\cosh$ function, we have:
$$ I = \frac{1}{\sqrt{2}}\sum_{n\geq 0}\int_{0}^{+\infty}\frac{dz}{\sqrt{z}\left[(2n+1)^2 \pi^2+ 4z^2\right]}=\frac{1}{2\sqrt{2\pi}}\sum_{n\geq 0}\frac{1}{(2n+1)^{3/2}}$$
hence:
$$ \boxed{I = \color{red}{\frac{2\sqrt{2}-1}{8\sqrt{\pi}}\,\zeta\left(\frac{3}{2}\right)}\approx 0.336859119428877} $$
For the two-variables integral you may ust switch to polar coordinates and apply the same technique, leading to a value of $\color{red}{\frac{\pi}{8}}$.
A: \begin{align}
\left(1 + a\mathrm{e}^{x^{2}}\right)^{-1}
&= \frac{1}{a}\mathrm{e}^{-x^{2}} \frac{1}{1+\frac{1}{a\mathrm{exp}(x^{2})}} \\
&= \frac{1}{a}\mathrm{e}^{-x^{2}} \sum\limits_{n=0}^{\infty} (-1)^{n} \left( \frac{1}{a}\mathrm{e}^{-x^{2}}\right)^{n} \\
&= \sum\limits_{n=0}^{\infty} (-1)^{n} \frac{1}{a^{n+1}} \mathrm{e}^{-(n+1)x^{2}} \\
\end{align}
Now integrate the exponential function
\begin{equation}
\int\limits_{0}^{\infty} \mathrm{e}^{-(n+1)x^{2}} \mathrm{d}x = \frac{\sqrt{\pi}}{2} \frac{1}{\sqrt{n+1}}
\end{equation}
Now we have
\begin{equation}
\int\limits_{0}^{\infty} \left(1 + a\mathrm{e}^{x^{2}}\right)^{-1} \mathrm{d}x =
\frac{\sqrt{\pi}}{2} \sum\limits_{n=1}^{\infty} (-1)^{n-1} \frac{1}{n^{1/2}} \frac{1}{a^{n}}
\end{equation}
Taking the derivative of both sides with respect to $a$:
\begin{equation}
\int\limits_{0}^{\infty} \mathrm{e}^{x^{2}} \left(1 + a\mathrm{e}^{x^{2}}\right)^{-2} \mathrm{d}x =
\frac{\sqrt{\pi}}{2} \sum\limits_{n=0}^{\infty} (-1)^{n-1} \frac{1}{n^{-1/2}} \frac{1}{a^{n+1}}
\end{equation}
Taking $\lim_{a \to 1}$
\begin{align}
\int\limits_{0}^{\infty} \mathrm{e}^{x^{2}} \left(1 + \mathrm{e}^{x^{2}}\right)^{-2} \mathrm{d}x &=
\frac{\sqrt{\pi}}{2} \sum\limits_{n=1}^{\infty} (-1)^{n-1} \frac{1}{n^{-1/2}} \\
&= \frac{\sqrt{\pi}}{2} \eta(-1/2) \\
&\approx 0.336859119
\end{align}
$\eta(s)$ is the Dirichlet eta function.
Addendum: Response to robjohn's comment.
From the reference I provided for the Dirichlet eta function:
"The following relation holds:
    $$\eta(s) = (1 − 2^{1 − s})\zeta(s)$$
While the Dirichlet series expansion for the eta function is convergent only for any complex number $s$ with real part > 0, it is Abel summable for any complex number. This serves to define the eta function as an entire function."
A: $$
\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}
$$
