built-in function to invert the matrix M , so the evaluation of the right-hand-side of Eq. (24) is easy. This procedure is encapsulated in a subroutine we called solvelinsys() . 

Before going on to the program for the iterative loop we need to decide when the gaps are small enough to stop the iteration. Thus we defined a test function which returns a boolean True if any of the four gaps, Δ g and Δ f as well as their slopes Δ g $^{′}$, and Δ f $^{′}$, is greater than δ . Again, in pseudocode, 

subroutine testgaps() boolean p1 = Abs[ggap] > delta boolean p2 = Abs[fgap] > delta boolean p3 = Abs[gpgap] > delta boolean p4 = Abs[fpgap] > delta boolean q = p1 OR p2 OR p3 OR p4 return q end subroutine 

The iteration will stop when q = False . 

The coding for running the iteration from an initial choice of parameters ˜ E and a 0 is as follows: 

{newtildeE,newa0} = {start_tildeE, start_a0} print starting parameters {newtildeE, newa0, a1} iter = 0 q = True do loop iter = iter+1 {oldtildeE,olda0} = {newtildeE,newa0} shootinandout(oldtildeE,olda0) gapsandpartials(oldtildeE,olda0) solvelinsys() // returns improved values of newtildeE and newa0 print iter and gaps print newtildeE and newa0 if q continue to next iteration, i.e., go to the top of the do loop else stop and break out of the do loop end if end do loop 

For the starting parameters used above, ˜ E = 0 . 82 and a 0 = 0 . 2, we find that our program closes the gaps in four iterations. The final parameter values are ˜ E = 0 . 727102 and a 0 = 0 . 194709. That ˜ E correponds, in more conventional units, to an energy of E = 0 . 306 GeV. Figure 2 shows plots of the final, converged g ( x ) and f ( x ) after normalization. Note the relatively large size of the lower component f ( x ), showing the importance of relativity in this case of massless quarks. 

