C FIGURE 1.6.1
      SUBROUTINE RKF(FSUB,N,T0,TFINAL,HINIT,U,TOL)                      
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
C                                                                       
C     ARGUMENT DESCRIPTIONS                                             
C                                                                       
C     FSUB   - (INPUT) THE NAME OF A USER-SUPPLIED SUBROUTINE TO EVALUATE
C              THE RIGHT HAND SIDE(S) OF THE DIFFERENTIAL EQUATIONS 
C              SYSTEM U' = F(T,U).  FSUB SHOULD HAVE THE CALLING SEQUENCE      
C                       SUBROUTINE FSUB(N,T,U,F)                        
C              AND SHOULD EVALUATE F(1)...F(N), GIVEN N,T AND U(1)...U(N).
C              FSUB MUST BE LISTED IN AN EXTERNAL STATEMENT IN THE      
C              CALLING PROGRAM.                                         
C     N      - (INPUT) THE NUMBER OF DIFFERENTIAL EQUATIONS
C     T0     - (INPUT) THE INITIAL VALUE OF T.                          
C     TFINAL - (INPUT) THE FINAL VALUE OF T.                            
C     HINIT  - (INPUT) THE INITIAL STEPSIZE TO BE USED.                 
C     U      - (INPUT) A VECTOR OF LENGTH N, HOLDING THE INITIAL VALUES 
C              OF U AT T=T0.                                            
C              (OUTPUT) THE SOLUTION VALUES AT T=TFINAL.                
C              THE SOLUTION COMPONENTS ARE PRINTED EACH TIME STEP.      
C     TOL    - (INPUT) THE DESIRED BOUND ON THE LOCAL ERROR PER UNIT TIME
C              OF THE SOLUTION COMPONENTS.                              
C                                                                       
      DIMENSION U(N),V1(N),V2(N),V3(N),V4(N),V5(N),V6(N),U4(N),U5(N),
     & UTEMP(N)
      TK = T0                                                           
      H = HINIT                                                         
    5 TKP1 = MIN(TK+H , TFINAL)                                         
      H = TKP1-TK                                                       
C                              TAKE A STEP WITH STEPSIZE H              
      CALL FSUB(N , TK , U , V1)                                        
      DO 10 I=1,N                                                       
   10 UTEMP(I) = U(I) + H*V1(I)/4                                       
      CALL FSUB(N , TK+H/4 , UTEMP , V2)                                
      DO 15 I=1,N                                                       
   15 UTEMP(I) = U(I) + H*(3*V1(I)+9*V2(I)) / 32                        
      CALL FSUB(N , TK+3*H/8 , UTEMP , V3)                              
      DO 20 I=1,N                                                       
   20 UTEMP(I) = U(I) + H*(1932*V1(I)-7200*V2(I)+7296*V3(I)) / 2197     
      CALL FSUB(N , TK+12*H/13 , UTEMP , V4)                            
      DO 25 I=1,N                                                       
   25 UTEMP(I) = U(I) + H*(8341*V1(I)-32832*V2(I)+29440*V3(I)           
     &                     -845*V4(I)) / 4104                           
      CALL FSUB(N , TK+H , UTEMP , V5)                                  
      DO 30 I=1,N                                                       
   30 UTEMP(I) = U(I) + H*(-6080*V1(I)+41040*V2(I)-28352*V3(I)          
     &                     +9295*V4(I)-5643*V5(I)) / 20520              
      CALL FSUB(N , TK+H/2 , UTEMP , V6)                                
      E = 0.0                                                           
      DO 35 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            
   35 E = MAX(E , ABS(U4(I)-U5(I))/H)                                   
C                              IF ERROR TOO LARGE, REPEAT STEP WITH H/2 
      IF (E .GT. TOL) THEN                                              
         H = H/2                                                        
         GO TO 5                                                        
      ENDIF                                                             
C                              OTHERWISE, UPDATE TK AND UK              
      TK = TKP1                                                         
      DO 40 I=1,N                                                       
         U(I) = U5(I)                                                      
   40 CONTINUE
      PRINT 45, TK,(U(I),I=1,N)                                         
   45 FORMAT (1X,5E15.5,/,(16X,4E15.5))                                 
C                              IF ERROR MUCH SMALLER THAN REQUIRED,     
C                              DOUBLE VALUE OF H FOR NEXT STEP.         
      IF (E .LT. 0.01*TOL) H = 2*H                                      
      IF (TK.LT.TFINAL) GO TO 5                                         
      RETURN                                                            
      END