It is best to vectorize your function using element-wise operations, then create it as a function of the two variables, then use integral2 to numerically integrate it:
m = 3;
n = 5;
p = 7;
f = @(theta,phi) m*p*cos(theta).^3.*sin(phi).^2 - n*p*cos(phi).^2.*cos(theta).^3 - m*n*cos(theta).*sin(phi).^2.*sin(theta).^2 + m*n*cos(phi).*cos(theta).*sin(phi).*sin(theta).^2;
f_int = integral2(f, 0, 2*pi, 0, 2*pi)
I created values for the constants to test the code.
