If you want to integrate a function of 4 variables, f(theta1,theta2,w,z), try iterating it as a single integral of a triple integral like so.
Q1 = @(theta1)integral3(@(theta2,w,z)f(theta1,theta2,w,z),theta2min,theta2max,wmin,wmax,zmin,zmax);
Now you can polish it off by
intgrl = integral(Q1,theta1min,theta1max,'ArrayValued',true);
The 'ArrayValued',true part ensures that integral() will not try to evaluate its integrand function, Q1, at anything other than scalar values of theta1. I'm not saying that this is the fastest way, but the 'ArrayValued' option is the easiest way because otherwise you have to do the work to define your integrand function, f, in such a way that it satisfies the requirements of integral or quadgk, namely that if given arrays of the same size as inputs for all its arguments, it returns an array of the same size. But quadgk returns a scalar, not an array, so you would need to write loops or use arrayfun to define your integrand function, and it can get complicated, or at least complicated-looking.