I want to evaluate the following integral numerically:
$$\int^{\infty}_{0}d\lambda_1\int^{\lambda_1}_{0}d\lambda_2 \lambda_1\lambda_2(\lambda_1-\lambda_2)\mathrm{Exp}[-\lambda_1\lambda_2-\lambda_1^2-\lambda_2^2]$$
The limits of integration impose the constraint that $\lambda_1 \geq \lambda_2$. The integral is something that can be solved analytically, and the result is
$$\frac{1}{18}(9-5\sqrt{3})\sqrt{\pi} = 0.0334547$$
I can't numerically evaluate this integral directly since there's a limit of $\lambda_1$ in one of the integrals, and symbols + numerical integration don't mix. So I change variables using $x_1=\lambda_1$, $x_2=\lambda_1-\lambda_2$. The Jacobian of this transform is $-1$, and the new integral is:
$$\int^\infty_0 dx_1\int^\infty_0 dx_2 x_1x_2(x_1-x_2)\mathrm{Abs}[-1]\mathrm{Exp}[-x_1(x_1-x_2)-x_1^2-(x_1-x_2)^2]$$
My argument for using these limits of integration are that $x_1=\lambda_1$ is a trivial transform, so it has the same limits of integration. Meanwhile since $\lambda_1 \geq \lambda_2$, their difference is positive, but could still range from $0$ to $\infty$.
The problem is, when I numerically integrate this with Mathematica, it produces the answer $-1.739$. I don't see how this is possible. Every term in the integrand is positive, so the integral cannot evaluate to a negative number. Even more bafflingly, if I integrate $x_2$ from $0$ to $-\infty$ (as opposed to $0$ to $\infty$), Mathematica produces the correct answer of $0.0334547$!
Can anyone explain what's going on? My best guess is that the Jacobian being negative flips the orientation of the integral such that I'm supposed to integrate from $0$ to $-\infty$, but that still doesn't explain how the limits of integration can change so haphazardly. After all, integrating from $0$ to $\infty$ is in general not related to integrating from $0$ to $-\infty$.
Your limits for $x_2$ look wrong. They should be from $x_1 \ to\ 0$, which need to be reversed.