4d integral in matlab

I'm trying to evaluate the integral using MATLAB with an equation containing 4 related random variables, thus the boundaries of the integrals are not constant.

There are 2 exponential pdfs, and 2 other hyper exponential , and 1 Reyleigh CDF all multiplied together and with (x - y - z).

I'm trying to evaluate it using integral Q = (@(w) integral3(@(x,y,z,w),xmin,xmax,ymin,ymax,zmin,zmax),wmin,wmax);

I'm getting an error always. this is my code down here :

u_x = 10; %  rate!
x_th = .3;
sigma = 1.33;
u_y = 10;
u_w = 100;
a= 1;
fun = @(x,y,z,w) (x - y - z )*u_x*exp(-u_x*x)*u_y*exp(-u_y*y)*((a/(a+1))*(a*u_w)*exp(-a*u_w*w)+((1/(a+1))*(u_w/a))*exp(-u_w*w/a))*((a/(a+1))*(a*u_w)*exp(-a*u_w*z)+((1/(a+1))*(u_w/a))*exp(-u_w*z/a))*(1-exp(-x_th/sigma^2))

xmin = @(y)y;
xmax = @(y,w)y + w;
ymin = 0;
ymax = inf;
zmin = 0;
zmax = @(w) w;
wmin = 0;
wmax = inf;

Q = integral(@(w) integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax),wmin,wmax);

ERROR MESSAGE :

Error using integral3 (line 63) XMIN must be a floating point scalar.

Error in numerical_int>@(w)integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax)

Error in integralCalc/iterateScalarValued (line 314) fx = FUN(t);

Error in integralCalc/vadapt (line 132) [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 83) [q,errbnd] = vadapt(@AToInfInvTransform,interval);

Error in integral (line 88) Q = integralCalc(fun,a,b,opstruct);

Error in numerical_int (line 28) Q = integral(@(w) integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax),wmin,wmax);

Solution

integral3() requires all the lower bounds to be real number
int() does not, but need to use syms function instead of function handle
For element-wise calculation use the dot. operator followed by the concerned operator like

1)     * ---> .*           2)     / ---> ./                3) ^ ---> .^

Read this for more information on how to use int() for nd integral

The code is as follows

syms x y z w
u_x = 10; %  rate!
x_th = .3;
sigma = 1.33;
u_y = 10;
u_w = 100;
a= 1;

fun = (x - y - z ).*u_x.*exp(-u_x*x).*u_y.*exp(-u_y.*y)...
    .*((a./(a+1)).*(a.*u_w).*exp(-a.*u_w.*w)+((1./(a+1)).*(u_w./a))...
    .*exp(-u_w.*w./a)).*((a./(a+1)).*(a*u_w)*exp(-a.*u_w.*z)...
    +((1./(a+1)).*(u_w./a))*exp(-u_w.*z./a)).*(1-exp(-x_th./sigma.^2));

xmin = y;
xmax = y + w;
ymin = 0;
ymax = inf;
zmin = 0;
zmax =  w;
wmin = 0;
wmax = inf;

% Integrate along x
intx = int(fun, x, xmin, xmax);

% Integrate along y
intxy = int(intx, y, ymin, ymax);

% Integrate along z
intxyz = int(intxy, z, zmin, zmax);

% Integrate along w
intxyzw = int(intxyz, w, wmin, wmax);

value = vpa(intxyzw, 3);

% 2.14e-5