In case anyone later comes along to read this post, I think I figured out a way to do the double integral using trapz.
% integration limits
x = linspace(-2,2);
y = linspace (0,3);
% given function
f = (x.^4) - (3*x.*y) + (6.*(y.^2));
q = x.^4 - x;
% implement trapezoidal rule
I = trapz(x,q.*trapz(y,f))
I don't know if it is right, but hopefully it will get you on the right track.