[theta, phi] = meshgrid(-pi:0.1:pi);
S=subplot(3,3);
S(1,1)=surf(ylm(1,0, -pi:0.1:pi, -pi:0.1:pi));
S(1,2)=surf(ylm(1,1, -pi:0.1:pi, -pi:0.1:pi));
S(1,3)=surf(ylm(2,1, -pi:0.1:pi, -pi:0.1:pi));
S(2,1)=surf(ylm(2,0, -pi:0.1:pi, -pi:0.1:pi));
S(2,2)=surf(ylm(2,1, -pi:0.1:pi, -pi:0.1:pi));
S(2,3)=surf(ylm(2,2, -pi:0.1:pi, -pi:0.1:pi));
S(3,1)=surf(ylm(3,0, -pi:0.1:pi, -pi:0.1:pi));
S(3,2)=surf(ylm(3,1, -pi:0.1:pi, -pi:0.1:pi));
S(3,3)=surf(ylm(3,2, -pi:0.1:pi, -pi:0.1:pi));