The case of massless quarks exhibits the role of relativity in a maximal way, meaning that the lower component f ( x ) is comparable in size to the upper component g ( x ). We will, later, choose an appropriate mass m for the massive quark flavors, s , c , and b , and solve for those wave functions. In these cases, as the quark mass m increases, f ( x ) becomes smaller relative to g ( x ), reflecting the transition to a non-relativistic limit where f ( x ) vanishes. Also, as will be seen, with increasing mass m , the heavier quark radial wave functions are progressively narrower than in the massless case. Thus we will, in Sec. IV, also include the effects of a color Coulomb-like vector potential V$_{v}$ coming from one-gluon exchange. This should prevent a too-rapid fall-off of the quark-quark interaction energies as the quark masses increase [17], but that is a matter for a separate paper. It simplifies the equations if we convert the ODEs to dimensionless form, dividing through by κ and defining a dimensionless distance x = κr : where ˜ E = E/κ , V$_{s}$ ( x ) = x − x 0 , and ˜ m = m/κ are now also dimensionless. Numerical integrations of these equations can be done by one or another of the Runge- Kutta methods.[18] To integrate these two first-order ODEs one needs to specify two bound- ary conditions (BCs) at the place where one begins the integration. However, as the radial wave functions eventually need to satisfy a normalization condition, this will determine one of the BCs, for example, the scale of the wave functions at asymp- totically large distances. In our first attempt to solve these equations we thought to integrate them outward from the origin with the hope of adjusting the energy eigenvalue ˜ E to assure that both g ( x ) and f ( x ) fall off to zero asymptotically. The first problem in doing this is the singularity in the equations at x = 0. This singularity, however, is easily avoided by starting instead at a small value of x = ϵ away from the origin. 5 A. Shooting Outwards g $^{′}$( x ) + k x g ( x ) − ( ˜ E + V$_{s}$ ( x ) + ˜ m ) f ( x ) = 0 , (10) f $^{′}$( x ) − k x f ( x ) + ( ˜ E − V$_{s}$ ( x ) − ˜ m ) g ( x ) = 0 , (11) ∫ ∞ 0 x $^{2}$dx [ ψ 2 $_{a}$( x ) + ψ 2 $_{b}$( x )] = ∫ ∞ 0 dx [ g $^{2}$( x ) + f $^{2}$( x )] = 1 , (12)