Try this —
syms phi r theta r0 r1 phi0 phi1 theta0 theta1
sympref('AbbreviateOutput',false);
func(r,phi,theta) = r.^2.*sin(phi).*(r.^2.*cos(phi).^2 + r.*cos(theta).*sin(phi) + r.*sin(phi).*sin(theta))
intr = int(func, r, r0, r1)
intphi = int(intr, phi, phi0, phi1)
inttheta = int(intphi, theta, theta0, theta1)
int3 = simplify(inttheta, 500)
Make appropriate changes to get the desired result.
.