Trying to approximate $E[f(X)]$ - Woflram Alpha gives $E[f(X)] \approx \frac{1}{\sqrt{3}}$ but I get $E[f(X)] \approx 0$? Let $X \sim \mathcal{N}(\mu_X,\sigma_X^2) = \mathcal{N}(0,1)$. Let $f(x) = e^{-x^2}$. I want to approximate $E[f(X)]$.
Wolfram Alpha gives
\begin{align}
E[f(X)] \approx \frac{1}{\sqrt{3}}.
\end{align}
Using a Taylor expansion approach, and noting that $f''(0) = -2$, I get
\begin{align}
E[f(X)] &\approx f(\mu_X) + \frac{f''(\mu_X)}{2} \sigma_X^2 \\
& = f(0) + \frac{f''(0)}{2} \\
& = 1 + \frac{-2}{2} \\
& = 0.
\end{align}
Why does my approximation fail to match the Wolfram Alpha result? What can be done to fix it?
 A: There's no need to "approximate" when you can derive the exact value of $\mathbb{E}[f(X)]$ . Let us apply the Law of the Unconscious Statistician (LoTUS) to obtain :
\begin{align*}
\mathbb{E}[f(X)] &= \int_{-\infty}^{+\infty} e^{-x^2} \cdot \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{x^2}{2}\right)~dx\\
&= 2\int_0^{+\infty} \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{3x^2}{2}\right)~dx\\
&= 2\int_0^{+\infty} \frac{1}{\sqrt{2\pi}} \cdot\frac{1}{\sqrt{6z}} e^{-z}~dz\\
&= \frac{2}{2\sqrt{3\pi}} \int_0^{+\infty} e^{-z} z^{\frac{1}{2} -1} ~dz\\
&=\frac{1}{\sqrt{3\pi}}\cdot \Gamma\left(\frac{1}{2}\right)\\
&= \frac{1}{\sqrt{3\pi}}\cdot \sqrt{\pi}\\
&= \frac{1}{\sqrt{3}}
\end{align*}
Hope this helps. :)
A: You does not need an approximation here. Use properties of moment generating functions, $X$ is standard normal so $X^2$ is chisquared with one df, with moment generating function $M_{X^2}(t)=\frac1{\sqrt{1-2t}}$ (for $t<1/2$.) Then note that
$$\DeclareMathOperator{\E}{\mathbb{E}}M_X(t)=\E e^{t X}
$$ is the definition, so that
$$\E e^{-X^2}=M_{X²}(-1)=\frac1{\sqrt{1-2\cdot (-1)}}=\frac1{\sqrt{3}}
$$
We can check that in R with a fast simulation (always a good idea to do a simulation check):
 mean( exp(-rnorm(1E6)^2) )
[1] 0.5774847
 1/sqrt(3)
[1] 0.5773503

Answer in comments:

What about if  was not standard normal but normal with mean $_$
and variance $^2_$. Can your approach still be used or is it
specific to the case of a standard normal distribution?

It can still be used. I will not give full details. First, the easy case $X \sim \mathcal{N}(0,\sigma^2)$. Then $X=(\sigma Z)^2$ with $Z$ standard normal, so in the above argument you get the argument $-\sigma^2$ in place of $-1$ for the mgf (moment generating function.)  For the fully general case, see for instance Moment-generating function (MGF) of non-central chi-squared distribution  and work from there.
