I will assume here that your function f(x,y,z) is written so that x, y, and z can be any size (vectors, matrices, N-D, whatever) and f(x,y,z) will be computed "element-wise". If this isn't the case, you can use arrayfun to build it up
fnew = @(x,y,z)arrayfun(f,x,y,z);
You should try your function on some array inputs if this is not clear, e.g. evaluate f([1,2,3],[4,5,6],[7,8,9]). The result should match [f(1,4,7),f(2,5,8),f(3,6,9)].
Now define a function of y and z where y and z are assumed to be scalars :
fyz_scalar = @(y,z)integral(@(x)f(x,y*ones(size(x)),z*ones(size(x)),xmin,xmax);
Now "vectorize" the function so that y and z can be arrays of any size (as long as they are the same size).
fyz = @(y,z)arrayfun(fyz_scalar,y,z);