As noted in the Introduction, the analytic solutions of the Dirac equation for the hydrogen atom were found long ago, [8] and these solutions do give the fine-structure splitting as a relativistic effect. Nonetheless, it is an interesting numerical exercise to see if the program for the linear scalar potential outlined above (or something like it) can be applied to the hydrogen atom case, and with enough accuracy. We want to solve Eqs. (3) and (4) with V$_{v}$ ( r ) = − Zα/r . We will set Z = 1, as these solutions are somewhat more delicate. Heavier, hydrogen-like atoms with Z > 1 can be done in a similar manner. Converting to dimensionless equations by dividing through by m$_{e}$ (instead of κ as before), where now x = m$_{e}$r and ˜ E = E/m$_{e}$ . The energy E is better written as E = m$_{e}$ − B , where B is the binding energy of the level of interest (which has principal quantum number n and orbital angular momentum quantum number l ). Note that, because of the smallness of the binding energies, E is only slightly less than m$_{e}$ . Thus ˜ E is less than (but very close to) 1, which is why we factored out a minus sign from the third term in Eq. (36), relative to that in Eq. (4). We will use the binding energy ˜ B = B/m$_{e}$ instead of ˜ E as one of the two parameters to be determined in the matching of the inwards and outwards integrations. For asymptotically large x the ODEs reduce to The solutions of these asymptotic equations at x = x max are which will be used as the BCs (starting values) for the inwards integration with a 1 as another parameter to be determined later by the normalization condition. The coefficient 15 g $^{′}$( x ) + k x g ( x ) − (1 + ˜ E + Zα/x ) f ( x ) = 0 , (35) f $^{′}$( x ) − k x f ( x ) − (1 − ˜ E − Zα/x ) g ( x ) = 0 , (36) g ( x max ) = a 1 e − µx max , f ( x max ) = − a 1 ( 1 − ˜ E 1 + ˜ $_{E}$) 1 / 2 e − µx max , (38) A. Boundary Conditions for Inward Integration g $^{′}$( x ) = (1 + ˜ E ) f ( x ) , f $^{′}$( x ) = (1 − ˜ E ) g ( x ) . (37)