The easiest way to integrate f(x,y,z,w) from a <= x <= b, c <= y <= d, e <= z <= g, and h <= w <= i is
Q = integral(@(x)integral3(@(y,z,w)f(x,y,z,w),c,d,e,g,h,i),a,b,'ArrayValued',true);
if the limits are functions of the preceding variables, i.e.:
c(x), d(x), e(x,y), g(x,y), h(x,y,z), i(x,y,z)
you will have to fix x in them as we did for the integrand:
Q = integral(@(x)integral3(@(y,z,w)f(x,y,z,w),c(x),d(x),@(y)e(x,y),@(y)g(x,y),@(y,z)h(x,y,z),@(y,z)i(x,y,z)),a,b,'ArrayValued',true);
The other, and perhaps faster way in some cases, would be to do an integral2(integral2(...)...) nesting, but this requires a bit of work to set up because of the vectorization that integral2 requires of the integrand.
