Return to Answer

Edit2

In Mathematica version 10, the volumetric slices can also be achieved using the ClipPlanes option:

With[{n = 3}, 
 Manipulate[
  Show[orb, Boxed -> False, Axes -> False, 
   ViewVector -> {{35, -35, 5}, {0, 0, 0}}, ViewAngle -> Pi/4, 
   ClipPlanes -> {{0, 1, 0, 1.2 y}}], {{y, 0}, -n^2, n^2}]]

The output is the same as in the previous animation.

Clear[rMin, rMax, r, θ, ϕ];
{rMin[n_, l_], rMax[n_, l_]} = 
  r /. Simplify[
    Solve[(l (l + 1))/r^2 - 2/r == -(1/n^2), r], 
    n > 0];

sphericalToCartesian = 
  Thread[{r, θ, ϕ} -> {Sqrt[x^2 + y^2 + z^2], 
     ArcCos[z/Sqrt[x^2 + y^2 + z^2]], Arg[x + I y]}];

ψ[n_, l_, m_][r_, θ_, ϕ_] := 
 Sqrt[ (n - l - 1)!/(n + l)!] E^(-(r/n)) ((2 r)/n)^l 2/
  n^2 LaguerreL[n - l - 1, 2 l + 1, (2 r)/n] SphericalHarmonicY[l, 
   m, θ, ϕ]

ClearAll[plotOrbital];
plotOrbital[f_, range_, contour_, opt : OptionsPattern[]] := 
 RegionPlot3D[
  Evaluate[Abs[f[r, \[Theta], \[Phi]] /. sphericalToCartesian]^2 > 
    contour], {x, -range, range}, {y, -range, range}, {z, -range, 
   range}, opt, Mesh -> False, PlotPoints -> 35, PlotStyle -> Orange, 
  Lighting -> "Neutral"]

grid = Table[
   Labeled[plotOrbital[ψ[n, l, m], 2 n^2, .05/n^6], 
    Row[{"n = ", n, ", ℓ = ", l, ", m = ", m}]], {n, 1, 
    3}, {l, 0, n - 1}, {m, 0, l}];

TableForm[grid]

Clear[rMin, rMax, r, θ, ϕ];
{rMin[n_, l_], rMax[n_, l_]} = 
  r /. Simplify[
    Solve[(l (l + 1))/r^2 - 2/r == -(1/n^2), r], 
    n > 0];

sphericalToCartesian = 
  Thread[{r, θ, ϕ} -> {Sqrt[x^2 + y^2 + z^2], 
     ArcCos[z/Sqrt[x^2 + y^2 + z^2]], Arg[x + I y]}];

ψ[n_, l_, m_][r_, θ_, ϕ_] := 
 Sqrt[ (n - l - 1)!/(n + l)!] E^(-(r/n)) ((2 r)/n)^l 2/
  n^2 LaguerreL[n - l - 1, 2 l + 1, (2 r)/n] SphericalHarmonicY[l, 
   m, θ, ϕ]

ClearAll[plotOrbital];
plotOrbital[f_, range_, contour_, opt : OptionsPattern[]] := 
 RegionPlot3D[
  Evaluate[Abs[f[r, θ, ϕ] /. sphericalToCartesian]^2 > 
    contour], {x, -range, range}, {y, -range, range}, {z, -range, 
   range}, opt, Mesh -> False, PlotPoints -> 35, PlotStyle -> Orange, 
  Lighting -> "Neutral", PlotTheme -> "Classic"]

grid = Table[
   Labeled[plotOrbital[ψ[n, l, m], 2 n^2, .05/n^6], 
    Row[{"n = ", n, ", ℓ = ", l, ", m = ", m}]], {n, 1, 
    3}, {l, 0, n - 1}, {m, 0, l}];

TableForm[grid]

Clear[rMin, rMax, r, θ, ϕ];
{rMin[n_, l_], rMax[n_, l_]} = 
  r /. Simplify[
    Solve[(l (l + 1))/r^2 - 2/r == -(1/n^2), r], 
    n > 0];

sphericalToCartesian = 
  Thread[{r, θ, ϕ} -> {Sqrt[x^2 + y^2 + z^2], 
     ArcCos[z/Sqrt[x^2 + y^2 + z^2]], Arg[x + I y]}];

ψ[n_, l_, m_][r_, θ_, ϕ_] := 
 Sqrt[ (n - l - 1)!/(n + l)!] E^(-(r/n)) ((2 r)/n)^l 2/
  n^2 LaguerreL[n - l - 1, 2 l + 1, (2 r)/n] SphericalHarmonicY[l, 
   m, θ, ϕ]

plotOrbital[f_, range_, contour_] := 
 RegionPlot3D[
  Evaluate[Abs[f[r, θ, ϕ] /. sphericalToCartesian]^2 > 
    contour], {x, -range, range}, {y, -range, range}, {z, -range, 
   range}, Mesh -> False, PlotPoints -> 35, PlotStyle -> Orange, 
  Lighting -> "Neutral"]

grid = Table[
   Labeled[plotOrbital[ψ[n, l, m], 2 n^2, .05/n^6], 
    Row[{"n = ", n, ", ℓ = ", l, ", m = ", m}]], {n, 1, 
    3}, {l, 0, n - 1}, {m, 0, l}];

TableForm[grid]

plotOrbital does the plotting for a given wave function f. In producing the grid of orbitals, I choose the parameters as described above.

Clear[rMin, rMax, r, θ, ϕ];
{rMin[n_, l_], rMax[n_, l_]} = 
  r /. Simplify[
    Solve[(l (l + 1))/r^2 - 2/r == -(1/n^2), r], 
    n > 0];

sphericalToCartesian = 
  Thread[{r, θ, ϕ} -> {Sqrt[x^2 + y^2 + z^2], 
     ArcCos[z/Sqrt[x^2 + y^2 + z^2]], Arg[x + I y]}];

ψ[n_, l_, m_][r_, θ_, ϕ_] := 
 Sqrt[ (n - l - 1)!/(n + l)!] E^(-(r/n)) ((2 r)/n)^l 2/
  n^2 LaguerreL[n - l - 1, 2 l + 1, (2 r)/n] SphericalHarmonicY[l, 
   m, θ, ϕ]

ClearAll[plotOrbital];
plotOrbital[f_, range_, contour_, opt : OptionsPattern[]] := 
 RegionPlot3D[
  Evaluate[Abs[f[r, \[Theta], \[Phi]] /. sphericalToCartesian]^2 > 
    contour], {x, -range, range}, {y, -range, range}, {z, -range, 
   range}, opt, Mesh -> False, PlotPoints -> 35, PlotStyle -> Orange, 
  Lighting -> "Neutral"]

grid = Table[
   Labeled[plotOrbital[ψ[n, l, m], 2 n^2, .05/n^6], 
    Row[{"n = ", n, ", ℓ = ", l, ", m = ", m}]], {n, 1, 
    3}, {l, 0, n - 1}, {m, 0, l}];

TableForm[grid]

plotOrbital does the plotting for a given wave function f. In producing the grid of orbitals, I choose the parameters as described above.

Edit: volumetric slices

To show that the RegionPlot approach is also able to provide multiple contours, and to show how these layered contours can give us information about the volume in 3D, here I first combine several outputs of plotOrbital with different contour values in a single Show. The surfaces have been given their own PlotStyle. Then I put the resulting 3D graphics in a Manipulate that allows you to slice through the layers in real time:

With[{n = 3, l = 2, m = 0}, 
 orb = Show[
   Table[plotOrbital[ψ[n, l, m], 2 n^2, c/n^6, 
     PlotStyle -> Directive[Opacity[0.6`], Hue[c]]], {c, 0.3`, 
     0.05`, -0.04`}], PlotRange -> All, BoxRatios -> Automatic]]

With[{n = 3},
 Manipulate[
  Show[orb, Boxed -> False, Axes -> False, 
   ViewVector -> {{35, -35, 5}, {0, 0, 0}}, ViewAngle -> Pi/4, 
   PlotRange -> {{-#, #}, {y, #}, {-#, #}} &[2 n^2]], {y, -2 n^2, 
   2 n^2}]]

The trick to get slices and hollow spaces between layers in RegionPlot3D is to cut off the plot with a reduced PlotRange. To still keep the object fixed in the view port (instead of moving around to stay centered with the changing PlotRange), I add fixed values for the ViewVector and ViewAngle. This trick using PlotRange doesn't work if you cut the range off directly inside the RegionPlot3D. You have to do it after the fact with Show.

Again, the drawing of multiple layers is very fast using RegionPlot3D. It could even be sped up more using ParallelTable.