For starting values (BCs) we can take advantage of the expected power-law behaviors of ψ$_{a}$ ( x ) and ψ$_{b}$ ( x ) near the origin for their known orbital angular momenta l and l $^{′}$. For example, if j = l + 1 $_{2}$, for which k = − ( l + 1) and the lower component’s l ′ = l + 1, one can choose as motivated by the non-relativistic limit, where ψ$_{a}$ → Ψ NR ∼ r l at small r . Substituting this f ( ϵ ) into the second ODE allows us to determine b 0 : The parameter a 0 at this point is arbitrary for now and will be fixed later by the normal- ization condition, Eq. (12). The problem with simply shooting outwards is what to choose for the energy E . The ODE solver for almost all initial guesses for E will have the well-known problem that the calculated g ( r ) and f ( r ) very soon blow up exponentially, either positively or negatively. It might be possible to iteratively refine the guess for E to eventually find decaying solution, but we soon decided to try a different method. A better way of assuring decaying solutions is to start from an asymptotic distance, where that behavior is built-in, and integrate the ODEs inward toward the origin. At large distances the ODEs simplify, for this potential, to which have solutions at large x = x max , Note that this asymptotic behavior is the same as that of the non-relativistic simple harmonic oscillator wave functions.[16] These forms will be used as starting values (BCs) for the inward integration of the full ODEs from x max . Again, a 1 can be taken as arbitrary for now, to 6 g ( ϵ ) = a 0 ϵ l +1 , f ( ϵ ) = b 0 ϵ l +2 , (13) f ( ϵ ) = − a 0 ( ˜ E − V$_{s}$ ( ϵ ) − ˜ m ) ϵ l $^{+2}$/ (2 l + 3) . (14) g $^{′}$( x ) − xf ( x ) = 0 , f $^{′}$( x ) − xg ( x ) = 0 , (15) g ( x max ) = a 1 e − x 2 $_{max}$/ 2 , f ( x max ) = − a 1 e − x 2 $_{max}$/ 2 . (16) B. Shooting Inwards