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, 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: 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. 11 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 {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