<< Graphics`ParametricPlot3D`

dzeta = 
      1. 
r = (Sqrt[(x^2) + (y^2)
 ])
c1s = Exp[-dzeta*r]
N1s = (Sqrt[dzeta^(3/(3  π))])
N2s = (Sqrt[dzeta^(3/(32 π))])
N2p = dzeta*N2s
N3d = (Sqrt[((dzeta^7*2)/(81^2* π))])
f1s = N1s*c1s
c2s = Exp[-dzeta*r/2]
f2s = N2s *(2 - dzeta*r)*c2s
f2px = N2p*x*c2s
c3d = Exp[-dzeta*r/3]
f3dxy = N3d*x*y*c3d
f3dx2 = 1/2*N3d*(x^2 - y^2)*c3d
c3dz = Exp[-dzeta*Sqrt[x^2 + z^2]/3]
f3dz2 = 1/(2*Sqrt[3])*N3d*(2*z^2 - x^2)*c3dz
ContourPlot[f1s, {x, -5, 5}, {y, -5, 5}, 
      PlotPoints -> 61]
Plot3D[f1s, {x, -2, 2}, {y, -2, 2}, 
      PlotPoints -> 61]
DensityPlot[f1s, {x, -2, 2}, {y, -2, 2}, 
      PlotPoints -> 61]
SphericalPlot3D[
      Abs[SphericalHarmonicY[0, 0, theta, phi]], {theta, 0, Pi, Pi/30}, {phi, 
        0, 2  Pi, Pi/35}, Boxed -> False, 
      Axes -> false] 
ContourPlot[f2s, {x, -10, 10}, {y, -10, 10}, 
      PlotPoints -> 61]
Plot3D[f2s, {x, -10, 10}, {y, -10, 10}, 
      PlotPoints -> 61]
DensityPlot[f2s, {x, -10, 10}, {y, -10, 10}, 
      PlotPoints -> 61]
SphericalPlot3D[
      2*Abs[SphericalHarmonicY[0, 0, theta, phi]], {theta, 0, Pi, 
        Pi/30}, {phi, 0, 2  Pi, Pi/15}, Boxed -> False, 
      Axes -> false] 
ContourPlot[f2px, {x, -10, 10}, {y, -10, 10}, 
      PlotPoints -> 61]
Plot3D[f2px, {x, -10, 10}, {y, -10, 10}, 
      PlotPoints -> 61]
DensityPlot[f2px, {x, -10, 10}, {y, -10, 10}, 
      PlotPoints -> 61]
SphericalPlot3D[
      Abs[SphericalHarmonicY[1, 0, theta, phi]], {theta, 0, Pi, Pi/30}, {phi, 
        0, 2  Pi, Pi/15}, Boxed -> False, 
      Axes -> false] 
ContourPlot[f3dxy, {x, -20, 20}, {y, -20, 20}, 
      PlotPoints -> 61]
Plot3D[f3dxy, {x, -20, 20}, {y, -20, 20}, 
      PlotPoints -> 61]
DensityPlot[f3dxy, {x, -20, 20}, {y, -20, 20}, 
      PlotPoints -> 61]
fxy = 
      Sqrt[15/4/Pi]*Sin[theta]*Cos[theta]*
        Cos[phi]
fxy1 = fxy/Cos[phi]*Sin[phi]
SphericalPlot3D[fxy, {theta, 0, Pi, Pi/30}, {phi, 0, 2  Pi, Pi/15}, 
      Boxed -> False, Axes -> false] 
SphericalPlot3D[fxy1, {theta, 0, Pi, Pi/30}, {phi, 0, 2  Pi, Pi/15}, 
      Boxed -> False, Axes -> false] 
ContourPlot[f3dx2, {x, -20, 20}, {y, -20, 20}, 
      PlotPoints -> 61]
Plot3D[f3dx2, {x, -20, 20}, {y, -20, 20}, 
      PlotPoints -> 61]
DensityPlot[f3dx2, {x, -20, 20}, {y, -20, 20}, 
      PlotPoints -> 61]
fxx = 
      Sqrt[15/16/Pi]*(Sin[theta])^2*
        Cos[2*phi]
SphericalPlot3D[fxx, {theta, 0, Pi, Pi/30}, {phi, 0, 2  Pi, Pi/15}, 
      Boxed -> False, Axes -> false] 
fxx1 = fxx/Cos[2*phi]*Sin[2*phi]
SphericalPlot3D[fxx1, {theta, 0, Pi, Pi/30}, {phi, 0, 2  Pi, Pi/15}, 
      Boxed -> False, Axes -> false] 
SphericalPlot3D[fxx1, {theta, 0, Pi, Pi/60}, {phi, 0, 2  Pi, Pi/15}, 
      Boxed -> False, Axes -> false]
SphericalPlot3D[fxx1, {theta, 0, Pi, Pi/60}, {phi, 0, 2  Pi, Pi/35}, 
      Boxed -> False, Axes -> false]
ContourPlot[f3dz2, {x, -25, 25}, {z, -25, 25}, 
      PlotPoints -> 61]
Plot3D[f3dz2, {x, -25, 25}, {z, -25, 25}, 
      PlotPoints -> 61]
DensityPlot[f3dz2, {x, -25, 25}, {z, -25, 25}, 
      PlotPoints -> 61]
SphericalPlot3D[
      Abs[SphericalHarmonicY[2, 0, theta, phi]], {theta, 0, Pi, Pi/30}, {phi, 
        0, 2  Pi, Pi/15}, Boxed -> False, 
      Axes -> false]