For the 2 s 1 2 state, starting with initial parameters B = 3 . 2 eV (presumably low), a 0 = 0 . 0007, and (fixed) a 1 = − 10 . 0, the program converges in four iterations to give final B = 3 . 40085 eV and a 0 = 0 . 00071905. The 2 s wave functions are plotted in Figs. 11 and 12. As expected on general grounds, unlike the ground state, here ψ$_{a}$ ( x ) and ψ$_{b}$ ( x ) each have one node. For the 2 p 1 2 state, starting with initial parameters B = 3 . 2 eV (also presumably low), a 0 = 0 . 0004, and (fixed) a 1 = 100 . 0, the program converged in four iterations to give final B = 3 . 40086 eV and a 0 = 0 . 000393546. In fact, the 2 s 1 2 and 2 p 1 2 levels are degenerate [27] as a result of an O (4) symmetry hiding in the equations. The small difference in the converged B ’s we find here is within the numerical precision of our Mathematica program. However, the 2 s 1 2 and 2 p 1 2 wave functions are very different . The 2 p 1 2 wave functions are shown in Figs. 13 and 14. Note that, here also, ψ$_{a}$ and ψ$_{b}$ each have one node. For the 2 p 3 2 state, starting this time with initial parameters B = 3 . 2 eV, a 0 = 0 . 000015, and a 1 = 100 . 0, the program converged in four iterations to give final B = 3 . 40144 eV and a 0 = 0 . 0000141876. The difference we find between the j = 3 2 and j = 1 2 energy levels, 0.58 meV, is the spin-orbit splitting, an intrinsically relativistic effect. The analytic value of this splitting is 0.453 meV [9, 27] and the difference here is due to the limited machine precision used in our Mathematica program. The 2 p 3 2 wave functions are shown in Figs. 15 and 16. The ψ$_{a}$ ( x )’s for the 2 p 3 2 and 2 p 1 2 states are very similar , but the lower component ψ$_{b}$ ( x )’s are quite different, reflecting the fact the lower components and the difference are both entirely due to relativity. For the 2 p 3 2 state also, ψ$_{a}$ and ψ$_{b}$ each have one node. There are no big surprises here, as Table I shows. From Eq. (39), the exponential fall-off of the wave functions is faster, as µ ∼ √$_{B}$ ∼ Z , so we adjust the values of x match and x max accordingly by dividing the hydrogenic values by Z . Running the code for the 1 s ground state for V$_{v}$ ( r ) = Zα/r works well for Z ≤ 100. As Z grows larger, the relativistic corrections to the energy become increasingly important. As Z approaches 1 /α = 137 . 036, however, the numerical solutions become inaccurate and, eventually, unstable. This can be seen in the rapid increase in the value of the parameter a 0 . The value of Z = 1 /α is where the Klein paradox [28] comes into play, since the analytic 18 E. Hydrogen-like Atoms with Z > 1