ClearAll[detail, R, Y, \[Psi], orb];
detail[n_, l_, m_] := 8 Total@{n, l, m}
F[n_, l_, x_ ] := 
  x^l Exp[-x/2] Factorial[n + 1] LaguerreL[n - l - 1, 2 l + 1, x];
R[n_ , l_, r_, 
   a_] := (2/(n a))^2 Sqrt[ a/4 (n - l - 1)!/((n + 1)!)^3] F[n, l, (
    2 r)/(n a)];
Y = SphericalHarmonicY;
\[Psi][n_ , l_, m_, r_, a_, \[Theta]_, \[Phi]_] := 
 R[n, l, r, a] Abs@Y[l, m, \[Theta], \[Phi]]
orb[n_, l_, m_, x_, y_, z_] := Module[{r = Norm@{x, y, z}},
   \[Psi][n, l, m, r, 1(*a*), ArcCos[z/r], ArcTan[x, y]]];
orbDensity[n_, l_, m_] := DensityPlot[orb[n, l, m, x, 0 (* doing a 2d cut in y=0 *), z],
  {x, -detail[n, l, m], detail[n, l, m]}, {z, -detail[n, l, m], 
   detail[n, l, m]}, 
  Mesh -> False, Frame -> False, PlotPoints -> 45, 
  ColorFunctionScaling -> True, ColorFunction -> "SunsetColors"]
MatrixForm@
 Table[MatrixForm@
   Table[Row@Table[{{n, l, m}, orbDensity[n, l, m]}, {m, 0, l}], {l, 
     0, n - 1}], {n, 1, 3}]