Maple3dplot

restart;
with(plots);
Y21 := (15/8)*cos(theta)^2*sin(theta)^2/Pi;
plot(Y21, theta = 0 .. Pi);
plot3d(Y21, phi = 0 .. 2*Pi, theta = 0 .. Pi, coords = spherical, scaling = constrained, grid = [100, 100]);
y00 := 1/(4*Pi);
plot3d(y00, phi = 0 .. 2*Pi, theta = 0 .. Pi, coords = spherical, scaling = constrained);
y10 := 3*cos(theta)^2/(4*Pi);
plot3d(y10, phi = 0 .. 2*Pi, theta = 0 .. Pi, coords = spherical, scaling = constrained);
y11 := 3*sin(theta)^2/(8*Pi);
plot3d(y11, phi = 0 .. 2*Pi, theta = 0 .. Pi, coords = spherical, scaling = constrained);