Small Increments Make the System Linear

My first attempt to calculate the time history of a geochemical system (Section 2.3) used the obvious approach (the direct Euler method) of evaluating the time derivatives and stepping forward. But it was not sue- 

In the expression for yp, substitute y x + delx) = y(jc) + dely to derive a system of algebraic equations for dely. Here are the equations presented earlier  

To clarify the notation I will write the indexes outside the parentheses thus 

Then this system of simultaneous linear algebraic equations can be solved using the subroutine GAUSS developed in Section 3.3. Because I have dropped the nonlinear term, I must always use a delx sufficiently small to ensure that all the dely values are indeed much smaller than the y values. 

Program DGC05 implements this solution with a subroutine to evaluate the coefficients of the sleq matrix. It calls subroutine GAUSS to calculate the increments dely, and then it steps forward in time and repeats. 

---