How to rewrite this integral for numerical evaluation

Question

$$\int_0^{\infty}\frac{e^{-z}}{\sqrt{z}}*e^{-\int_\alpha^{\infty}(v-\int_0^{\infty}\frac{vue^{-\frac{u^2}{2}}}{1+2uvz}du)dv}dz$$

I'm trying to numerically evaluate this integral in Matlab using quad and dblguad functions, but I'm running into problems due to the way the inner integrals are structured. I cant use dblquad (double integration) to integrate over $u$ and $v$ because the function in the innermost integral is a function of $u$, $v$ and $z$. I guess my question is there any way I can get $z$ out of the two inner integrals. This way I can evaluate the inner integrals and all is left is just a function of $z$ for the outermost integral.

This is my code, but it's not working, because of the problem I mentioned above.

F1 = @(u,v) u.*v.*exp(-0.5.*u.^2)./(1+2*z.*u.*v);  
F2 = @(v) v;
F3 = @(z) exp(-z)./sqrt(z);
I1 = dblquad(F1,0,1e5,2,1e5);
I2 = quad(F2,2,1e5);
quad(F3*exp(-(I2-I1)),0,1e5);

Answer 1

First observe that

\begin{align} \int_a^\infty \int_0^\infty v e^{-u} - \frac{v u e^{-u^2/2}}{1 + 2 u v z} \, du \, dv &= \int_a^\infty v \int_0^\infty e^{-u} \, du - \int_0^\infty \frac{v u e^{-u^2/2}}{1 + 2 u v z} \, du \, dv \\ &= \int_a^\infty v \cdot 1 - \int_0^\infty \frac{v u e^{-u^2/2}}{1 + 2 u v z} \, du \, dv. \end{align}

This allows us to rewrite the inner two integrals as a standard double integral:

$$\int_0^{\infty}\frac{e^{-z}}{\sqrt{z}} \exp\left[-\int_a^\infty \left( v - \int_0^\infty \frac{vue^{-\frac{u^2}{2}}}{1+2uvz}\,du \right) \,dv\right]\,dz = \int_0^{\infty}\frac{e^{-z}}{\sqrt{z}} \exp\left[-\int_a^\infty \int_0^\infty v e^{-u} - \frac{vue^{-\frac{u^2}{2}}}{1+2uvz}\,du\,dv\right]\,dz.$$

Next, in MATLAB we define a function myfun.m which, given a vector of $z$-values, evaluates the function

$$z \mapsto \frac{e^{-z}}{\sqrt{z}} \exp\left[-\int_a^\infty \int_0^\infty v e^{-u} - \frac{vue^{-\frac{u^2}{2}}}{1+2uvz}\,du\,dv\right].$$

function f = myfun(z)

a = 1;
max_x = 1000;
max_y = max_x;

f = zeros(size(z));

for idx = 1:numel(z)
    f(idx) = integral2(@(u, v) v.*exp(-u) - v.*u.*exp(-u.^2/2)./(1 ...
                               + 2*u.*v*z(idx)),0,max_x,a,max_y);
end

f = exp(-z - f)./sqrt(z);

The reason to use max_x and max_y instead of inf is because otherwise the integral diverges and MATLAB throws an error. (This might be an indication that the problem is not well-posed; I'm not sure.) It might be possible to vectorize the loop, but it isn't clear to me how to do that.

Lastly, q = integral(@myfun, 0, inf) will compute the desired integral (1).