This turned out to be surprisingly tricky, but here is the solution. First, define a function
function fzv = integrateOverSecondPair(fun,x,y,zmin,zmax,vmin,vmax)
for ii=numel(x):-1:1
g = @(z,v) fun(x(ii),y,z,v);
fzv(ii) = dblquad(g,zmin,zmax,vmin,vmax);
end
Then use the commands
fzv = @(x,y) integrateOverSecondPair(@f,x,y,zmin,zmax,vmin,vmax);
fxyzv = dblquad(fzv,xmin,xmax,ymin,ymax);
Best Answer