plo[n_, l_, m_] := DensityPlot[Evaluate[Abs[psi[n, l, m, r, u, v]]^2 //. 
       {r -> Sqrt[x^2 + z^2], u -> ArcCos[z/r], v -> 0}], 
     {x, -10*(1.5*n + 1), (Plus[10])*(1.5*n + 1)}, {z, -10*(1.5*n + 1), 
      (Plus[10])*(1.5*n + 1)}, ColorFunction -> Hue, Mesh -> False, 
     PlotPoints -> 200, PlotLabel -> {n, l, m}]
 
psi[n_, l_, m_, r_, theta_, phi_] := radialFunction[n, l, r]*
     SphericalHarmonicY[l, m, theta, phi]
 
radialFunction[n_, l_, r_] := (2*Sqrt[(n - l - 1)!/(n + l)!^3]*
       f[n, l, (2*r)/(n*a)])/n^2/a^(3/2)
 
f[n_, l_, x_] := x^l*Exp[-(x/2)]*(n + l)!*LaguerreL[n - l - 1, 2*l + 1, x]
 
a = 1
 
disp := Display["t.ps", #1, "EPS"] & 
