It turns out that this iterative procedure works very well, with the gaps in the slopes of g ( x ) and f ( x ) at the match point also going to zero. At this point we should remark that others might have done this calculation but by zeroing the gaps in the logarithmic derivatives, such as( dg ( x ) /dx ) /g ( x ). This has the advantage of removing the scale dependence of the functions g and f , but our choice of defining the gaps as relative, as in Eq. (18), is essentially equivalent, as it also involves (semi-)local differences divided by the (semi-)local value. We implemented the iterative numerical process discussed above by developing Mathe- matica notebooks, [19] but any reasonable programming language could be used instead. An example notebook is available from our group’s web site.[20] (We especially invite our Canadian colleagues to convert this to a Maple program!) Here, as we describe the calcu- lation of the ground state radial wave functions for the GMSS potential of Eq. (9), we will comment on a number of the programming issues that we encountered. The first thing to be done is to set the value of k in the equations to be solved. It is an integer and depends upon the angular momentum quantum numbers, Eq. (5). For the 1 s , j = 1 / 2 ground state, l = 0, l ′ = 1, and thus k = − 1. As noted earlier, we work in units where ¯ h = c = 1 [21] so that the dimensions of Eqs. (7) and (8) are fm − $^{1}$. For the cases discussed here and below, we have fixed the constants in the GMSS potential at κ 2 = 0 . 9 GeV / fm (i.e., κ = 2 . 14 fm − $^{1}$), and r 0 = 0 . 705 fm. Making the equations dimensionless, as mentioned after Eq. (11), this linear scalar potential simplifies to ˜ V$_{s}$ ( x ) = x − x 0 , where x 0 = κr 0 = 1 . 506. For dimensionless distances we chose ϵ = 10 − $^{6}$, x match = 1 . 0, and x max = 6 . 0. We also define a small increment δ = 0 . 0001, which will be used when we calculate the partial derivatives in Eq. (22) numerically and for testing when the gaps to be zeroed are small enough. The two integrations, inwards and outwards, were done using Mathematica’s NDSolve function, but one could use any standard Runge-Kutta procedure.[18] To proceed, we need the four boundary conditions (BCs) from Eq. (13) and Eq. (16). We have defined a subroutine, shootinandout , to do this. This subroutine is, of course, a function of the 9 E. The 1 s Radial Wave Functions and Some Programming Details