How to evaluate this integral in Matlab?

I'm trying to evaluate the following integral in matlab:

Here is my code:

alpha = 2;
F1 = @(u,v) 2*u.*v.*exp(-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);

I'm getting the errors shown below. These errors don't show much, but I'm guessing it's because of the way I wrote F1. I defined F1 as a function of u and v for the double integral, but there is also a variable z which is the variable of the outer integral. I did that because there is no way I can separate z from the inner integrals. Is there any better way to write this integration?

Undefined function or variable 'z'
Error in ==> @(u,v)2*u.*v.*exp(-u.^2)./(1+2*z.*u.*v)
Error in ==> dblquad>innerintegral at 73
fcl = intfcn(xmin, y(1), varargin{:}); %evaluate only to get the class below
Error in ==> quad at 76
y = f(x, varargin{:});
Error in ==> dblquad at 53
Q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ...

I'm choosing 1e5 to represent infinity.

Answer 1

[Mike Hosea]

Sorry I missed this question back in October. It is not easy to do these problems because one has to keep track of what works with array input and what requires scalar input. You also have to build the answer from the inside out. Here is how I would approach it:

upper = inf;
% Define the innermost integrand integral as a function of u, v and z.
F1 = @(u,v,z)2.*u.*v.*exp(-u.*u)./(1 + 2.*u.*v.*z);
% Define the inner integral as a function of v and z.  It is
% important to realize that this only works for scalar v and z.
i1scalar = @(v,z)integral(@(u)F1(u,v,z),0,upper);
% Define the middle integrand as a function of v and z.  Because i1scalar
% only works for scalar v and z, this function only works for scalar v and
% z.
F2scalar = @(v,z)v - i1scalar(v,z);
% Enhance the function to accept vector v (but not z).
F2 = @(v,z)arrayfun(@(v)F2scalar(v,z),v);
% Define the middle integral as a function of z and alpha, which must be
% scalars.
i2scalar = @(z,alpha)integral(@(v)F2(v,z),alpha,upper);
% Define the integrand for the outer integration as a function of z and
% alpha.  Both must be scalars
Fscalar = @(z,alpha)exp(-z)./sqrt(z).*exp(-i2scalar(z,alpha));
% Enhance the integrand to accept a vector z input (but only scalar alpha).
F3 = @(z,alpha)arrayfun(@(z)Fscalar(z,alpha),z);
% Define the final result as a function of alpha.
Qscalar = @(alpha)integral(@(z)F3(z,alpha),0,upper);
% Enhance Qscalar to accept vector alpha.
Q = @(alpha)arrayfun(Qscalar,alpha);

You can avoid the use of arrayfun, if you like, with the 'ArrayValued' option, but that is usually slower.

Unfortunately, the middle integral is a big problem. Here v ranges from alpha to infinity, so the way this integrand is written, relying on cancellation, we will lose all precision there. The problem will need to be reformulated in order to be integrated successfully. I suppose you can substitute something finite for the upper limit. I made it easy for you to try this, but this often results in spurious zeros from numerical integration routines, and I wouldn't trust the answer.