In this notebook, we will use Mathematica to solve the Schrödinger equation for hydrogen atom and plot some of the results. You will need the following two packages.

Off[General::spell]; Off[General::spell1];
Needs["Calculus`VectorAnalysis`"]
Needs["Graphics`ParametricPlot3D`"]

The time-independent Schrödinger equation is written as:
	- ((ℏ^2)/(2  μ))(∇^( 2))ψ((Overscript[r,  →]) ) + V((Overscript[r,  →]) ) ψ((Overscript[r,  →]) ) = E ψ((Overscript[r,  →]) )

Typically the potential is a function of r only, i.e., V((Overscript[r,  →]) )=V(r), it is natural to choose the spherical coordinates (r,θ,ϕ).

SetCoordinates[
        Spherical[r, θ, ϕ]];
CoordinateRanges[]

ClearAll["Global`*"]

Let's look for solutions that are separable into products:
	Ψ(r,θ,ϕ)= R(r) Θ(θ) Φ(ϕ)
For hydrogen, the potential is just the Coulomb's potential:
	V(r)= ((- (e^2))/(4  π (Subscript[ϵ, 0]) r)),
and the Schrödinger equation is given as follows:

ψ = 
        R[r] Θ[θ] Φ[ϕ];
eqn = -((ℏ^2)/(2  μ)) Laplacian[ψ] - 
        ((e^2)/(4 π (Subscript[ϵ, 0]) r)) ψ - ℰ ψ

The Schrödinger equation can be separated into three ordinary differential equations with two constants ℓ(ℓ+1) and (m^2):

eqn1 = Expand[eqn/ψ];
eqn1 = Expand[(r^2) eqn1];
rEqn = 
        eqn1[[1]] + 
          eqn1[[2]] + 
          eqn1[[3]] + 
          eqn1[[5]] + ℓ (ℓ + 1) ((ℏ^2)/(2  μ));
rEqn = 
        Simplify[
          rEqn ((2  μ)/(ℏ^2))];
rEqn = Expand[R[r] rEqn]
angEqn = 
        eqn1[[4]] + 
          eqn1[[6]] + 
          eqn1[[7]] - ℓ (ℓ + 1) ((ℏ^2)/(2  μ));
angEqn = 
        Expand[Sin[θ]^2 angEqn];
angEqn = 
        Expand[(-((2  μ)/(ℏ^2))) angEqn];
θEqn = 
        angEqn[[1]] + 
          angEqn[[2]] + 
          angEqn[[3]] + 
          angEqn[[4]] - 
          (m^2);
θEqn = 
      Expand[Θ[θ] θEqn]
ϕEqn = 
        angEqn[[5]] + 
          (m^2);
ϕEqn = Expand[Φ[ϕ] ϕEqn]

The ϕ-equation and its solution:

The solution for the ϕ-equation is just (e^(i m ϕ)). The boundary condition Φ(0) = Φ(2π) requires that 𝓂 be integers.

ϕAns = 
      DSolve[ϕEqn ==  0, Φ[ϕ], ϕ]
ϕAns = 
        TrigToExp[Φ[ϕ] /. ϕAns];
ϕAns = ϕAns /. (Sqrt[(m^2)]) ->  
            m;
ϕAns = 
        Collect[ϕAns, {(E^(-I m ϕ)), (E^(I m ϕ))}];
(Subscript[Φ, m])_[ϕ_] = ϕAns[[1, 
          1]] /. (C[1]/2 - 
              (1/2) I C[2]) ->  1

The θ-equation and its solution:

With z = cosθ, the θ-equation can be written in terms of z as:
	(1-(z^2))\(TraditionalForm\`\(\(d\^2\) \(Θ(z)\)\)\/\(d\ z\^2\)\)- 2 z ((d Θ(z))/(d z)) + [ℓ(ℓ+1) - ((m^2)/(1 - (z^2)))] Θ(z) = 0
which is just the associated Legendre equation.

zEqn = θEqn //. {Cos[θ] ->  
              z, 
                          Θ'[θ] ->  (∂_θ) 
                  Cos[θ] (∂_z) Θ[
                    z], 
                          Θ''[θ] ->  (∂_(θ, θ))Cos[θ] (∂_z) Θ[
                      z] + 
                                         ((((∂_θ) Cos[θ]))^2) (∂_(z, z))Θ[
                      z], 
                          Sin[θ]^2 ->  (1  - 
                  (z^2)), 
                          Θ[θ] ->  Θ[
                z]};
zEqn = 
      Simplify[Collect[
          zEqn, 
                                            {(∂_(z, z))Θ[z], (∂_z) Θ[z], Θ[
              z]}]]
zEqn = 
      zEqn /. (-1 + 
              (z^2)) ->  -Sin[z]^2
zEqn = 
      Collect[Expand[zEqn/Sin[z]^2], Θ[
          z]]
zEqn = 
      zEqn //. {Sin[z]^2 ->  (1 - 
                (z^2)), 
                           Csc[z]^2 ->  (1/(1 - (z^2)))}

The general solution for the associated Legendre equation is given in terms of the Legendre functions of the first and second kind. For integer 𝓂, the Legendre functions of the second kind have logarithmic singularities at z=±1. So we retain only the first kind for Θ(z).

zAns = 
      DSolve[zEqn ==  0, Θ[z], 
        z]
(Subscript[Θ, ℓ_, m_])[z_] = 
      zAns[[1, 1, 2, 1]] /. 
        C[1] ->  1

The spherical harmonic is then defined in terms of Θ(z) and Φ(ϕ) with appropriate normation.

(Subscript[Y, ℓ_, m_])[θ_, ϕ_] := 
                           (Sqrt[(((2  ℓ + 1))/(4  π)) (((ℓ - Abs[m])!)/((ℓ + Abs[m])!))]) (Subscript[Θ, ℓ, m])[
          Cos[θ]] (Subscript[Φ, m])[ϕ]

The r-equation and its solution:

The r-equation can be simplified if we change variables:
	R(r) = u(r) / r.

rEqn1 = 
      rEqn //. {R[r] ->  
            u[r]/r, 
                              (∂_r) 
              R[r] ->  (r (∂_r) u[r] - u[r])/
              (r^2), 
                              (∂_(r, r))R[
                r] ->  
                                          (∂_(r, r))u[
                    r]/r - (2/(r^2)) (∂_r) 
                  u[r]
                                         + 
              2  u[r]/(r^3)}
rEqn1 = Expand[rEqn1/r]

Let's again change the variables:
	r = ρ / κ
	(Subscript[ρ,  0])= \(TraditionalForm\`\(μ\ e\^2\)\/\(2  π\ κ\ ℏ\^2\ ϵ\_\(\(\ \)\(0\)\)\)\)

ϱEqn = 
      rEqn1 /. {r ->  ϱ/κ, 
                         u[r] ->  
            u[ϱ], 
                         (∂_r) 
              u[r] ->  κ (∂_ϱ) 
                u[ϱ], 
                         (∂_(r, r))u[
                r] ->  (κ^2) (∂_(ϱ, ϱ))u[ϱ]}
ϱEqn = 
      Expand[ϱEqn/(κ^2)]
ϱEqn = ϱEqn /. {-2 ℰ μ/(ℏ^2) ->  (κ^2), 
                        - (((e^2) μ u[ϱ])/(2 π κ ϱ (ℏ^2) (Subscript[ϵ, 0]))) ->  -(Subscript[ϱ, 0]) u[ϱ]/ϱ}

We can peel off the asymptotic behaviour of u(ρ) by introducing the new function v(ρ):
	u(ρ) = (ρ^(ℓ + 1))(e^(-ρ))v(ρ)

ϱEqn=ϱEqn/.{u[ϱ] ->  (ϱ^(ℓ + 1)) (E^(-ϱ)) v[ϱ],
                        (∂_(ϱ, ϱ))u[ϱ]-> (E^(-ϱ)) ℓ (1 + ℓ) (ϱ^(-1 + ℓ)) v[ϱ]-2 (E^(-ϱ)) (1 + ℓ) (ϱ^ℓ) v[ϱ]+(E^(-ϱ)) (ϱ^(1 + ℓ)) v[ϱ]+2 (E^(-ϱ)) (1 + ℓ) (ϱ^ℓ) v'[ϱ]-2 (E^(-ϱ)) (ϱ^(1 + ℓ)) v'[ϱ]+(E^(-ϱ)) (ϱ^(1 + ℓ)) v''[ϱ]}
ϱEqn = 
      Simplify[ϱEqn]
ϱEqn = 
      Collect[ϱEqn, {v[ϱ], v'[ϱ], v''[ϱ]}]
ϱEqn = 
      Expand[ϱEqn/((E^(-ϱ)) (ϱ^ℓ))]
ϱEqn = 
      Collect[ϱEqn, {v[ϱ], v'[ϱ], v''[ϱ]}]

We are looking for a polynomial solution to the equation for v(ϱ). To perserve the function v(ϱ), let's temperorily replace the function v(ϱ) with w(ϱ)  and expand w(ϱ) as:
	w(ϱ) = (∑_(j = 0)^N) (Subscript[a, j]) (ϱ^j)

wEqn=ϱEqn/.{v[ϱ] ->  w[ϱ],v'[ϱ]-> w'[ϱ],v''[ϱ]-> w''[ϱ]}
w[ϱ_] =  (∑_(j = 0)^N) 
                (Subscript[a, j]) (ϱ^j);
expansion = (wEqn /. N ->  5) + 
        O[ϱ]^6
ans = 
      LogicalExpand[
        expansion ==  0]
coeff = 
      Solve[ans, {(Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3]), (Subscript[a, 4]), (Subscript[a, 5])}]

ϱEqn = ϱEqn /. (Subscript[ϱ, 0]) ->  
          2  n

ϱAns = 
      DSolve[ϱEqn ==  0, 
        v[ϱ], ϱ]
(Subscript[v, n_, ℓ_])[ϱ_] = ϱAns[[1, 1, 2, 2]] /. 
        C[2] ->  1

(Subscript[f, n_, ℓ_])[ϱ_] := (E^(-ϱ)) (((2  ϱ))^ℓ) (n + ℓ)! (Subscript[v, n, ℓ])[ϱ]
(Subscript[R, n_, ℓ_])[
        r_] := (((a^(-3/2)) 2 (Sqrt[(((n - ℓ - 1)!)/(((n + ℓ)!)^3))]) (Subscript[f, n, ℓ])[((2 r)/(n a))])/(n^2))

(Subscript[R, 1, 0])[r/2]
(Subscript[R, 2, 0])[r/2]
(Subscript[R, 2, 1])[r/2]
(Subscript[R, 3, 0])[r/2]
(Subscript[R, 3, 1])[r/2]
(Subscript[R, 3, 2])[r/2]

Plot[
        Evaluate[
          Table[(((r (Subscript[R, i, 0])[r/2]))^2), {i, 1, 4}] /. a -> 1], {r, 
          0, 25}, PlotRange -> {0, 0.55}, 
        AxesLabel -> {"r(a)", "P(r) (1/a)"}, 
        PlotStyle -> {{RGBColor[0, 0, 1]}, {RGBColor[1, 0, 
                1]}, {RGBColor[0, 1, 0]}, {RGBColor[1, 0.5, 0]}}, 
        ImageSize -> 400, 
        Epilog -> {Text["1 s", {3, 0.5}], 
            Text["2 s", {6, 0.22}], Text["3 s", {15, 0.14}], 
            Text["4 s", {24, 0.1}]}];
Plot[
        Evaluate[
          Table[(((r (Subscript[R, i, 1])[r/2]))^2), {i, 2, 4}] /. a -> 1], {r, 
          0, 25}, PlotRange -> {0, 0.3}, 
        AxesLabel -> {"r(a)", "P(r) (1/a)"}, 
        PlotStyle -> {{RGBColor[1, 0, 1]}, {RGBColor[0, 1, 
                0]}, {RGBColor[1, 0.5, 0]}}, ImageSize -> 400, 
        Epilog -> {Text["2 p", {4, 0.22}], 
            Text["3 p", {12, 0.13}], 
            Text["4 p", {24, 0.1}]}];
Plot[
        Evaluate[
          Table[(((r (Subscript[R, i, 2])[r/2]))^2), {i, 3, 4}] /. a -> 1], {r, 
          0, 25}, PlotRange -> {0, 0.2}, 
        AxesLabel -> {"r(a)", "P(r) (1/a)"}, 
        PlotStyle -> {{RGBColor[0, 1, 0]}, {RGBColor[1, 0.5, 0]}}, 
        ImageSize -> 400, 
        Epilog -> {Text["3 d", {9, 0.125}], 
            Text["4 d", {21, 0.08}]}];

(Subscript[Ψ, n_, ℓ_, m_])[r_, θ_, ϕ_] := 
      (Subscript[R, n, ℓ])[
          r/2] (Subscript[Y, ℓ, m])[θ, ϕ]
\(\(Ψ_\^*\) := Ψ /. 
        Complex[u_, v_] -> \ 
          Complex[u, \(-v\)]\)
\(⟨ψ_ \[VerticalSeparator] 
          A_ \[VerticalSeparator] φ_⟩ := ∫\_0\%∞\(∫\_0\%π\(∫\_0\%\(2  π\)\(ψ\^*\)\ A\ φ\ r\^2\ Sin[θ] \[DifferentialD]ϕ \[DifferentialD]θ \[DifferentialD]r\)\)\)

The average values of the various powers of r are given as:
	⟨nℓm| r |nℓm⟩ = (1/2)[3 (n^2)- ℓ(ℓ+1)] a
	⟨nℓm| (r^2) |nℓm⟩ = (1/2)(n^2)[5 (n^2)+1 - 3ℓ(ℓ+1)] (a^2)
	⟨nℓm| (r^(-1)) |nℓm⟩ = (1/((n^2) a))
	⟨nℓm| (r^(-2)) |nℓm⟩ = (2/((n^3)(2  ℓ + 1) (a^2)))
In the cell below we evaluate the average values for n=1, ℓ=0, and m=0.

⟨(Subscript[Ψ, 1, 0, 0])[
          r, θ, ϕ] \[VerticalSeparator] 
        r \[VerticalSeparator] (Subscript[Ψ, 1, 0, 0])[
          r, θ, ϕ]⟩
⟨(Subscript[Ψ, 1, 0, 0])[
          r, θ, ϕ] \[VerticalSeparator] 
        (r^2) \[VerticalSeparator] (Subscript[Ψ, 1, 0, 0])[
          r, θ, ϕ]⟩
⟨(Subscript[Ψ, 1, 0, 0])[
          r, θ, ϕ] \[VerticalSeparator] 
        1/r \[VerticalSeparator] (Subscript[Ψ, 1, 0, 0])[
          r, θ, ϕ]⟩
⟨(Subscript[Ψ, 1, 0, 0])[
          r, θ, ϕ] \[VerticalSeparator] 
        1/(r^2) \[VerticalSeparator] (Subscript[Ψ, 1, 0, 0])[
          r, θ, ϕ]⟩

hydrogenSphericalPlot3D[n_, ℓ_, 
        m_] := 
  SphericalPlot3D[
        Evaluate[
          Conjugate[(Subscript[Ψ, n, ℓ, m])[
                  a, θ, ϕ]] (Subscript[Ψ, n, ℓ, m])[
                a, θ, ϕ] /. 
            a ->  1], 
{θ, 
          0, π}, {ϕ, -π, π}, Boxed ->  False, 
        DisplayFunction -> Identity, 
        PlotLabel -> {n, ℓ, m}, 
 
        PlotPoints ->  45, Axes ->  None]

Show[
        GraphicsArray[
{{hydrogenSphericalPlot3D[4, 0, 
                0]}, {hydrogenSphericalPlot3D[4, 1, 
                0], 
hydrogenSphericalPlot3D[4, 1, 
                1]}, {hydrogenSphericalPlot3D[4, 2, 
                0], 
hydrogenSphericalPlot3D[4, 2, 
                1], 
hydrogenSphericalPlot3D[4, 2, 
                2]}, 
{hydrogenSphericalPlot3D[4, 3, 
                0], 
hydrogenSphericalPlot3D[4, 3, 
                1], 
hydrogenSphericalPlot3D[4, 3, 
                2], 
hydrogenSphericalPlot3D[4, 3, 
                3]}}, 
GraphicsSpacing -> 
            0.1]];

Clear[x, y, z]
hydrogenPlot[n_, ℓ_, m_] := 
      DensityPlot[
        Evaluate[
          Conjugate[(Subscript[Ψ, n, ℓ, m])[
                  r, θ, ϕ]] (Subscript[Ψ, n, ℓ, m])[
                r, θ, ϕ] //. {a -> 1, 
              r -> (Sqrt[(x^2) + (z^2)]), θ -> 
                ArcCos[(z/r)], ϕ -> 
                0}], {x, -(8/2) (1.5 n + 1), 
          (8/2) (1.5 n + 1)}, {z, -(8/2) (1.5 n + 1), 
          (8/2) (1.5 n + 1)}, Mesh -> False, 
        PlotPoints -> 200, FrameLabel -> {"x", "z"}, 
        PlotLabel -> {n, ℓ, m}, 
        DisplayFunction -> Identity, FrameTicks -> None, 
        ImageSize -> 600]

Show[
        GraphicsArray[
{{hydrogenPlot[4, 0, 
                0]}, {hydrogenPlot[4, 1, 0], 
              hydrogenPlot[4, 1, 1]}, {hydrogenPlot[4, 2, 0], 
              hydrogenPlot[4, 2, 1], 
              hydrogenPlot[4, 2, 2]}, {hydrogenPlot[4, 3, 0], 
              hydrogenPlot[4, 3, 1], 
              hydrogenPlot[4, 3, 
                2]}}, 
GraphicsSpacing -> 
            0.03]];

Show[
        GraphicsArray[{{hydrogenPlot[1, 0, 0]}, {hydrogenPlot[2, 0, 0], 
              hydrogenPlot[2, 1, 0], 
              hydrogenPlot[2, 1, 1]}, {hydrogenPlot[3, 0, 0], 
              hydrogenPlot[3, 1, 0], 
              hydrogenPlot[3, 1, 1]}, {hydrogenPlot[3, 2, 0], 
              hydrogenPlot[3, 2, 1], 
              hydrogenPlot[3, 2, 2]}, {hydrogenPlot[4, 0, 0], 
              hydrogenPlot[4, 1, 0], 
              hydrogenPlot[4, 1, 1]}, {hydrogenPlot[4, 2, 0], 
              hydrogenPlot[4, 2, 1], hydrogenPlot[4, 2, 2]}}, 
          GraphicsSpacing -> 0.03]];