be later fixed by the normalization condition, but we need to make an initial guess for the energy eigenvalue ˜ E . The integration will go from x max back to x = ϵ , again to avoid the singularity at x = 0. The problem with this shooting inwards method is much the same as that for shooting outwards. Unless one somehow is able to guess the exact value of ˜ E , the solutions will blow up exponentially as one nears the origin, and usually well before that. Thus, after experiencing these well-known problems with ODEs with eigenvalues, we con- cluded that we had to combine the two methods, invoking a well-known numerical method, “shoot and match.” In this case we shoot both outwards from x = ϵ to a point in the middle, x match , and inwards from x max back to x match . The two values for g out ( x match ) and g in ( x match ) will generally differ for a given choice of ˜ E and another parameter, which we choose to be a 0 . That is likewise the case for f out ( x match ) and f in ( x match ). What we therefore need is to find a way to iteratively vary ˜ E and a 0 to reduce the two gaps, Δ g and Δ f , to zero. [The asymptotic parameter a 1 is chosen to be fixed, i.e., not to be varied, but of a size that makes the initial determination of the gaps at x match reasonably small; it will in the end be determined by the normalization condition, Eq. (12).] We need to define the gaps Δ g , etc., at the match point. Initially, we simply chose them in terms of their actual differences, and likewise for Δ f , Δ g $^{′}$, and Δ f $^{′}$. Indeed, this often works for the problem we originally set out to do, namely, to calculate wave functions for the (nearly) massless u and d quarks. However, as the quark mass ˜ m increases, the system becomes more and more non-relativistic. This means that the lower component wave function, f ( x ), becomes relatively small com- pared with g ( x ). Thus a better definition for the gaps would be to make them relative to their average value, e.g., and likewise for Δ f ( ˜ E, a 0 ). 7 Δ g ( ˜ E, a 0 ) = 2 [ g out ( x match ) − g in ( x match ) ] g out ( x match ) + g in ( x match ) , (18) Δ g ( ˜ E, a 0 ) = g out ( x match ) − g in ( x match ) , (17) C. Shooting In and Out and Matching