How can I evaluate a Double Integral with an additional variable dependent integral contained inside using dblquad?

I am trying to compute an integral of the following form:

doubleintegral(g(x,y)*[singleintegral(f(x,y,z)dz)]dxdy).

I am able to compute this integral using quadgk and dblquad if f is only dependent on two variables, like so:

 g = @(x,y) x.*y;
 f = @(z,y) z.*y;
 Int1 = @(y) quadgk(@(z) f(z,y),0,1);
 IntF = dblquad(@(x,y) g(x,y).*Int1(y),0,1,0,1);

However, if I try to extend f to 3 variables, I cannot get the same algorithm to work:

 g = @(x,y) x.*y;
 f = @(z,y,x) z.*y.*x;
 Int1 = @(y,x) quadgk(@(z) f(z,y,x),0,1);
 IntF = dblquad(@(x,y) g(x,y).*Int1(y,x),0,1,0,1);

I am not sure what I am missing to extend this same process for f it is of 3 variables, as the above code gives me a number of errors, including something about the matrix dimensions not agreeing.

Any help with computing this would be greatly appreciated!

Answer 1

[Mike Hosea]

I forget DBLQUAD's arcane requirements on the input sizes. If memory serves, the integrand has to work when one is a vector and the other is a scalar. Your definition of Int1 is only going to return a scalar if it works at all, so it can't possibly return a vector of the same size as whichever input is a vector. The newer INTEGRAL2 function (and QUAD2D) instead require that the integrand work when both inputs are arrays of arbitrary size (the same size). The output needs to be an array of the same size. This is easy to accomplish with ARRAYFUN.

g = @(x,y) x.*y;
f = @(z,y,x) z.*y.*x;
Int1 = @(y,x) quadgk(@(z) f(z,y,x),0,1); % Only scalar y and x work.
Int1v = @(y,x)arrayfun(Int1,y,x); % Vector y and x now work.
IntF = integral2(@(x,y) g(x,y).*Int1v(y,x),0,1,0,1);