% FIGURE 1.6.1
      function U = RKF(FSUB,N,T0,TFINAL,HINIT,U,TOL)
%
%     ARGUMENT DESCRIPTIONS
%
%     FSUB   - (INPUT) STRING CONTAINING THE NAME OF A USER-SUPPLIED
%              FUNCTION TO EVALUATE THE RIGHT HAND SIDE(S) OF THE
%              DIFFERENTIAL EQUATIONS SYSTEM U' = F(T,U).
%              F = FSUB(N,T,U) SHOULD EVALUATE F(1)...F(N), GIVEN
%              N,T AND U(1)...U(N).
%     N      - (INPUT) THE NUMBER OF DIFFERENTIAL EQUATIONS
%     T0     - (INPUT) THE INITIAL VALUE OF T.
%     TFINAL - (INPUT) THE FINAL VALUE OF T.
%     HINIT  - (INPUT) THE INITIAL STEPSIZE TO BE USED.
%     U      - (INPUT) A VECTOR OF LENGTH N, HOLDING THE INITIAL VALUES
%              OF U AT T=T0.
%              (OUTPUT) THE SOLUTION VALUES AT T=TFINAL.
%              THE SOLUTION COMPONENTS ARE PRINTED EACH TIME STEP.
%     TOL    - (INPUT) THE DESIRED BOUND ON THE LOCAL ERROR PER UNIT TIME
%              OF THE SOLUTION COMPONENTS.
%
      TK = T0;
      H = HINIT;
      while (TK < TFINAL)
         TKP1 = min(TK+H , TFINAL);
         H = TKP1-TK;
%                              TAKE A STEP WITH STEPSIZE H
         V1 = feval(FSUB,N,TK,U);
         for I=1:N
            UTEMP(I) = U(I) + H*V1(I)/4;
         end
         V2 = feval(FSUB,N,TK+H/4,UTEMP);
         for I=1:N
            UTEMP(I) = U(I) + H*(3*V1(I)+9*V2(I)) / 32;
         end
         V3 = feval(FSUB,N,TK+3*H/8,UTEMP);
         for I=1:N
            UTEMP(I) = U(I) + H*(1932*V1(I)-7200*V2(I)+7296*V3(I)) / 2197;
         end
         V4 = feval(FSUB,N,TK+12*H/13,UTEMP);
         for I=1:N
            UTEMP(I) = U(I) + H*(8341*V1(I)-32832*V2(I)+29440*V3(I) ...
                            -845*V4(I)) / 4104;
         end
         V5 = feval(FSUB,N,TK+H,UTEMP);
         for I=1:N
            UTEMP(I) = U(I) + H*(-6080*V1(I)+41040*V2(I)-28352*V3(I) ...
                           +9295*V4(I)-5643*V5(I)) / 20520;
         end
         V6 = feval(FSUB,N,TK+H/2,UTEMP);
         E = 0.0;
         for I=1:N
            U4(I) = U(I) + H*(2375*V1(I)+11264*V3(I)+10985*V4(I) ...
                          -4104*V5(I)) / 20520;
            U5(I) = U(I) + H*(33440*V1(I)+146432*V3(I)+142805*V4(I) ...
                          -50787*V5(I)+10260*V6(I)) / 282150;
            E = max(E , abs(U4(I)-U5(I))/H);
         end
%                              IF ERROR TOO LARGE, REPEAT STEP WITH H/2
         if (E > TOL)
            H = H/2;
            continue
         end
%                              OTHERWISE, UPDATE TK AND UK
         TK = TKP1
         for I=1:N
            U(I) = U4(I);
         end
         U
%                              IF ERROR MUCH SMALLER THAN REQUIRED,
%                              DOUBLE VALUE OF H FOR NEXT STEP.
         if (E < 0.01*TOL)
             H = 2*H;
         end
      end