I. INTRODUCTION Dirac formulated relativistic quantum mechanics in the late 1920’s.[1] Since field theory had not yet been developed, the relativistic aspect led to a number of confusions related to currents since only total charge is conserved. With the recognition of the appearance of antiparticles, and the development of field theory, analogous to the transition from the Canonical to the Grand Canonical ensemble in statistical mechanics, these issues were re- solved. Nonetheless, the solution of the equation itself, ignoring these deeper aspects, has proven valuable even in more modern contexts. For example, the MIT bag model of quark confinement,[2] and the theories that evolved from it (chiral bag, cloudy bag, etc.)[3] depend on solving the Dirac equation as a wave function for a particle in an effective potential. Through boundary conditions and other approximations, the Dirac equation has even been employed in this way in the study of nuclear structure.[4, 5] In particle physics, the quark model has been employed, with great success, to describe baryon and meson states and their structure. Usually this is done in a non -relativistic model, [6] where the potentials follow the patterns expected by a Foldy-Wouthuysen reduction[7] of a simpler, theoretically motivated potential in the Dirac equation. In our case, however, we wanted to see how these structures and solutions appear without the non-relativistic reduction approximations, by solving the Dirac equation itself for the simpler potential with fewer adjustable parameters. Going back somewhat in time, very soon after Dirac’s initial formulation, Darwin [8] found an analytic solution (!) for the radial wave functions for hydrogen-like atoms in terms of confluent hypergeometric functions. Since then, solutions of the Dirac equation have been important in atomic physics, even to recent times.[9, 10] For example, more complicated central potentials than Zα/r , the fourth component of a Lorentz four-vector, usually require finding a numerical solution. Another example is a muonic atom, in which a muon is in an s-state about a heavy nucleus having a realistic charge distribution.[11] Finding such numerical solutions involves solving coupled ordinary differential equations for the upper and lower components of the Dirac wave function, [10, 12] 2 ψ$_{jlm}$ ( r ) =    ψ$_{a}$ ( r ) Y l jm − iψ$_{b}$ ( r ) Y l ′ jm    = 1 r    ig ( r ) Y l jm − f ( r ) Y l ′ jm    , (1)