Psi[r_, z_, n_, l_, m_]=Exp[-Sqrt[(r^2 + z^2)]/2]*Sqrt[r^2 + z^2]^l*LaguerreL[n - l - 1, 2  l + 1, Sqrt[r^2 + z^2]]*SphericalHarmonicY[l, m, ArcTan[r/z], 0]

Table[
      ContourPlot[
        Log[Abs[Psi[r, z, n, l, m]]^2], {r, -5 n, 5  n}, {z, -5 n, 
          5  n}, PlotPoints -> 50, ContourLines -> False, 
        FrameLabel -> {"ρ/a", "z/a"}, 
        PlotLabel -> {"nlm=", n, l, m}], {n, 1, 3}, {l, 0, 
        n - 1}, {m, 0, l}]