I have an algorithm for which I need to compute some multidimensional integrals. Sometimes, Integrate works but some others Integrate just fails to compute the integral.
A simple example:
$$ \int^{\infty}_{-\infty}\int^{\infty}_{-\infty} -\frac{2 x^2 \left(e^{2 y+1}+1\right) (\text{erf}(y)+\text{erf}(y+1)) e^{-x^2-(y+1)^2}}{\pi } {\rm d}x {\rm d}y$$
Calling integrate,
Integrate[-((2 E^(-x^2 - (1 + y)^2) (1 + E^(1 + 2 y)) x^2 (Erf[y] + Erf[1 + y]))/\[Pi]), {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]
yields
$$ \int_{-\infty }^{\infty } -\frac{e^{-(y+1)^2} \left(e^{2 y+1}+1\right) (\text{erf}(y)+\text{erf}(y+1))}{\sqrt{\pi }} \, {\rm d}y $$
which is 0 (Interestingly, WolframAlpha indeed handles this integral, while my Mathematica does not). Now, Calling NIntegrate,
NIntegrate[-((2 E^(-x^2 - (1 + y)^2) (1 + E^(1 + 2 y)) x^2 (Erf[y] +
Erf[1 + y]))/\[Pi]), {x, -Infinity, Infinity}, {y, -Infinity, Infinity}, AccuracyGoal -> 14]
yields $-5.2\times 10^{-16}$ (After a bunch of warnings). I was wondering I there was any way of further calling NIntegrate whenever Integrate fails. I know that I can call
N[Integrate[...]]
but I have not figured out how to pass method options to this. I am looking for a generic solution that works for higher dimensional integrals. If it helps, my integrals all involve exponentials and error functions similar to that one.
Thanks in advance.