c  This is a modified version of DVERK that uses a defect based error and
c stepsize control. It produces a continous approximation to the solution
c on each step. The continuous solution can be accessed on each step
c by using option C(9) to return control before accepting a step. On
c return from DVERK the driver program can check to ensure that the
c current step will be accepted, IND=5, and if so the continuous
c solution can be evaluated at any point in the current step
c (xt between X and XTRIAL) using the routine INTRP. A derivative
c evaluation routine, DERIV, is also available for use at this point.
c
c This routine uses the same calling sequence and shares the same
c overview as the original code DVERK. For this reason the internal
c documentation for DVERK is relevant to this modified version as
c well. In particular it describes the options available and how
c they can be invoked.
c
c Note that W must be declared to be W(n,14) and the value of YTRIAL
c (on return after an attempted step) is stored in W(*,13).
c
C YTRIAL IS IN W(K,13)
      SUBROUTINE DERIV( N, X, Y, XOUT, YPOUT, HTRIAL, NW, W)            
C
      INTEGER N, NW, I                                                  
C     EXTERNAL FCN                                                      
      DOUBLE PRECISION X,Y(N),XOUT,YPOUT(N),W(NW,14),TAU,HTRIAL         
      DOUBLE PRECISION D0,D1,D2,D3,D4,D5,D6
      TAU = ( XOUT - X ) / HTRIAL                                       
      D0 = -40.D0/(17.D0*HTRIAL) * (TAU - 1.D0)                         
     +     * ((((20.D0*TAU - 40.4D0)*TAU + 23.1D0)*TAU                  
     +        - 1.7D0)*TAU - .85D0)                                     
     +     - (TAU - 1.D0)**2/(17.D0*HTRIAL)                             
     +     *(((1.6D3*TAU-2424.D0)*TAU+924.D0)*TAU-34.D0)                
      D1 = -2.D3/8721.D0 * ((3.D0*TAU - 4.D0)*TAU + 1.D0)               
     +     * ((( 15.1D0*TAU - 32.542D0)*TAU                             
     +          + 21.903D0)*TAU - 4.3605D0)                             
     +     - TAU*(TAU - 1.D0)**2/8721.D0                                
     +     *((90600.D0*TAU - 130168.D0)*TAU + 43806.D0)                 
      D2 = 40.D0/(17.D0*HTRIAL) * TAU                                   
     +     * (((( 20.D0*TAU - 80.4D0)*TAU + 123.9D0)*TAU                
     +        - 88.3D0)*TAU + 25.65D0)                                  
     +     + TAU**2/(17.D0*HTRIAL)                                      
     +     *(((1.6D3*TAU-4824.D0)*TAU+4956.D0)*TAU-1766.D0)             
      D3 = 4.D3/102.D0 * (3.D0*TAU - 2.D0)*TAU                          
     +     * ((( 1.9D0*TAU - 3.698D0)*TAU                               
     +         + 2.08D0)*TAU - .2565D0)                                 
     +     + (TAU - 1.D0)*TAU**2 / 51.D0                                
     +     * ((114.D2*TAU - 14792.D0)*TAU + 4160.D0)                    
      D4 = -8.D3/459.D0 * (TAU*(TAU - 1.D0)*(2.D0*TAU - 1.D0))          
     +     * (( 2.3D0*TAU - 4.136D0)*TAU + 14.877D0/8.D0)               
     +     -(TAU*(TAU-1.D0))**2/918.D0 * (368.D2*TAU-33088.D0)          
      D5 = 250.D0/918.D0 * TAU*(TAU - 1.D0)*(2.D0*TAU - 1.D0)           
     +     * ((800.D0*TAU - 800.D0)*TAU + 57.D0)                        
     +     + 125.D0/918.D0*(TAU*(TAU-1.D0))**2 * (1.6D3*TAU-8.D2)       
      D6 = -16.D3/8721.D0 * TAU*(TAU - 1.D0)*(2.D0*TAU - 1.D0)          
     +     * ((2.D2*TAU - 2.D2)*TAU + 27.D0)                            
     +     - 8.D3/8721.D0 * (TAU*(TAU - 1.D0))**2                       
     +     * (4.D2*TAU - 2.D2)                                          
C                                                                       
      DO 40 I = 1,N                                                     
         YPOUT(I) = D0*Y(I) + D1*W(I,1) + D2*W(I,13) + D3*W(I,9)        
     +             + D4*W(I,10) + D5*W(I,11) + D6*W(I,12)               
 40   CONTINUE                                                          
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE INTRP( N, X, Y, XOUT, YOUT, HTRIAL, NW, W)
C
      INTEGER N, NW, I
      DOUBLE PRECISION X,Y(N),XOUT,YOUT(N),W(NW,14),TAU,HTRIAL          
      DOUBLE PRECISION D0,D1,D2,D3,D4,D5,D6
      TAU = ( XOUT - X ) / HTRIAL                                       
      D0 = -20.D0/17.D0 * (TAU-1.D0)**2                                 
     +     * ((((20.D0*TAU - 40.4D0)*TAU + 23.1D0)*TAU                  
     +        - 1.7D0)*TAU - .85D0)                                     
      D1 = -2.D3/8721.D0 * TAU * (TAU-1.D0)**2                          
     +     * ((( 15.1D0*TAU - 32.542D0)*TAU                             
     +          + 21.903D0)*TAU - 4.3605D0) * HTRIAL                    
      D2 = 20.D0/17.D0 * TAU**2                                         
     +     * (((( 20.D0*TAU - 80.4D0)*TAU + 123.9D0)*TAU                
     +        - 88.3D0)*TAU + 25.65D0)                                  
      D3 = 4.D3/102.D0 * (TAU-1.D0) * TAU**2                            
     +     * ((( 1.9D0*TAU - 3.698D0)*TAU                               
     +         + 2.08D0)*TAU - .2565D0) * HTRIAL                        
      D4 = -4.D3/459.D0 * TAU**2 * (TAU-1.D0)**2                        
     +     * (( 2.3D0*TAU - 4.136D0)*TAU + 14.877D0/8.D0) * HTRIAL      
      D5 = 25.D3/459.D0 * TAU**2 * (TAU-1.D0)**2                        
     +     * (( 2.D0*TAU - 2.D0)*TAU + .1425D0) * HTRIAL                
      D6 = -8.D5/8721.D0 * TAU**2 * (TAU-1.D0)**2                       
     +     * ((2.D0*TAU - 2.D0)*TAU + .27D0) * HTRIAL                   
C                                                                       
      DO 40 I = 1,N                                                     
         YOUT(I) = D0*Y(I) + D1*W(I,1) + D2*W(I,13) + D3*W(I,9)         
     +             + D4*W(I,10) + D5*W(I,11) + D6*W(I,12)               
 40   CONTINUE                                                          
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
c      SUBROUTINE DEFECT(N,XOUT,DFCT,FCN,X,Y,C,NW,W)                     
c      DOUBLE PRECISION XOUT,DFCT(N),X,Y(N),C(30),W(NW,16),YPP(10),      
c     +                     YOUT(10),YP(10)                              
c      INTEGER N,I,NW                                                    
c      EXTERNAL FCN                                                      
C                                                                       
C     CALL INTRP(N,X,Y,XOUT,YOUT,C(18),NW,W)
C     CALL FCN(  XOUT,YOUT,YP)
C     CALL DERIV(N,X,Y,XOUT,YPP,C(18),NW,W)
C     DO 10 I = 1,N                                                     
C        DFCT(I) = YPP(I) - YP(I)                                       
c10    CONTINUE                                                          
c      RETURN                                                            
c      END                                                               
C                                                                       
      SUBROUTINE DVERK (N, FCN, X, Y, XEND, TOL, IND, C, NW, W)         
      INTEGER N, IND, NW, K, I                                          
      DOUBLE PRECISION X, Y(N), XEND, TOL, C(30), W(NW,14), TEMP        
      DOUBLE PRECISION A1, A3, A4, A5, A6, A7, A8, A0,                  
     +     BY0, BY1, BY5, BF0, BF1, BF5,                                
     +     CY0, CY1, CY5, CF0, CF1, CF5,                                
     +     D0, D1, D2, D3, D4, D5, D6                                   
      PARAMETER ( A1 = 49.D0/640.D0,     A3 = 4375.D0/11968.D0,        
     +            A4 = 23.D0/384.D0,     A5 = -33.D0/1955.D0,              
     +            A6 = 0.D0,             A7 = 125.D0/8832.D0,          
     +            A8 = -43.D0/1408.D0,   A0 = 1.D0/32.D0 )              
      PARAMETER ( BY0 = 128.D0/3125.D0, BY1 = 2592.D0/3125.D0,          
     +            BY5 = 81.D0/625.D0,    BF0 = 18.D0/3125.D0,           
     +            BF1 = -162.D0/3125.D0, BF5 = 162.D0/3125.D0 )         
      PARAMETER ( CY0 = 5427.D0/400000.D0, CY1 = 380133.D0/400000.D0,   
     +            CY5 = 361.D0/10000.D0,     CF0 = 1539.D0/800000.D0,   
     +            CF1 = -29241.D0/800000.D0, CF5 = 3249.D0/200000.D0 )  
C                                                                       
C ADDED EXTERNAL STATEMENT FOR FCN TO AVOID A WARNING MESSAGE.          
C      EXTERNAL FCN                                                      
C***********************************************************************
C                                                                       
C     ******************************************************************
C     * BEGIN INITIALIZATION, PARAMETER CHECKING, INTERRUPT RE-ENTRIES *
C     ******************************************************************
C                                                                       
C  ......ABORT IF IND OUT OF RANGE 1 TO 6                               
      IF (IND.LT.1 .OR. IND.GT.6) GO TO 500                             
C                                                                       
C CASES - INITIAL ENTRY, NORMAL RE-ENTRY, INTERRUPT RE-ENTRIES          
      GO TO (5, 5, 45, 1111, 2222, 2222), IND                           
C CASE 1 - INITIAL ENTRY (IND .EQ. 1 OR 2)                              
C  .........ABORT IF N.GT.NW OR TOL.LE.0                                
    5    IF (N.GT.NW .OR. TOL.LE.0.D0) GO TO 500                        
      IF (IND.EQ. 2) GO TO 15                                           
C       INITIAL ENTRY WITHOUT OPTIONS (IND .EQ. 1)                      
C       SET C(1) TO C(9) EQUAL TO 0                                     
        DO 10 K = 1, 9                                                  
           C(K) = 0.D0                                                  
   10   CONTINUE                                                        
        GO TO 35                                                        
   15  CONTINUE                                                         
C       INITIAL ENTRY WITH OPTIONS (IND .EQ. 2)                         
C       MAKE C(1) TO C(9) NON-NEGATIVE                                  
        DO 20 K = 1, 9                                                  
           C(K) = DABS(C(K))                                            
   20   CONTINUE                                                        
C       MAKE FLOOR VALUES NON-NEGATIVE IF THEY ARE TO BE USED           
        IF (C(1).NE.4.D0 .AND. C(1).NE.5.D0) GO TO 30                   
           DO 25 K = 1, N                                               
              C(K+30) = DABS(C(K+30))                                   
   25      CONTINUE                                                     
   30   CONTINUE                                                        
   35  CONTINUE                                                         
C    INITIALIZE RREB, DWARF, PREV XEND, FLAG, COUNTS                    
       C(10) = 2.D0**(-56)                                              
       C(11) = 1.D-35                                                   
C    SET PREVIOUS XEND INITIALLY TO INITIAL VALUE OF X                  
       C(20) = X                                                        
       DO 40 K = 21, 24                                                 
          C(K) = 0.D0                                                   
   40  CONTINUE                                                         
       GO TO 50                                                         
C CASE 2 - NORMAL RE-ENTRY (IND .EQ. 3)                                 
C  .........ABORT IF XEND REACHED, AND EITHER X CHANGED OR XEND NOT     
   45  IF (C(21).NE.0.D0 .AND.                                          
     +      (X.NE.C(20) .OR. XEND.EQ.C(20))) GO TO 500                  
C     RE-INITIALIZE FLAG                                                
       C(21) = 0.D0                                                     
       GO TO 50                                                         
C CASE 3 - RE-ENTRY FOLLOWING AN INTERRUPT (IND .EQ. 4 TO 6)            
C    TRANSFER CONTROL TO THE APPROPRIATE RE-ENTRY POINT..........       
C    THIS HAS ALREADY BEEN HANDLED BY THE COMPUTED GO TO        .       
C END CASES       V                                                     
   50 CONTINUE                                                          
C                                                                       
C     END INITIALIZATION, ETC.                                          
C                                                                       
C      *****************************************************************
C     * LOOP THROUGH THE FOLLOWING 4 STAGES, ONCE FOR EACH TRIAL  STEP *
C     * UNTIL THE OCCURRENCE OF ONE OF THE FOLLOWING       *            
C     *    (A) THE NORMAL RETURN (WITH IND .EQ. 3) ON REACHING XEND IN *
C     *        STAGE 4       *                                          
C     *    (B) AN ERROR RETURN (WITH IND .LT. 0) IN STAGE 1 OR STAGE 4 *
C     *    (C) AN INTERRUPT RETURN (WITH IND  .EQ.  4,5  OR  6),  IF *  
C     *        REQUESTED, IN STAGE 1 OR STAGE 4        *                
C     ******************************************************************
C                                                                       
99999 CONTINUE                                                          
C                                                                       
C ***************************************************************       
C * STAGE 1 - PREPARE - DO CALCULATIONS OF  HMIN,  HMAX,  ETC., *       
C * AND SOME PARAMETER  CHECKING,  AND  END  UPWITH  SUITABLE *         
C * VALUES OF HMAG, XTRIAL AND HTRIAL IN PREPARATION FOR TAKING *       
C * AN INTEGRATION STEP.        *                                       
C ***************************************************************       
C                                                                       
C***********ERROR RETURN (WITH IND=-1) IF NO OF FCN EVALS TOO GREAT     
       IF (C(7).EQ.0.D0 .OR. C(24).LT.C(7)) GO TO 100                   
          IND = -1                                                      
          RETURN                                                        
  100  CONTINUE                                                         
C                                                                       
C           STORE K1 IN THE FIRST COLLUMN OF W. UNLESS THIS IS THE FIRST
C           STEP THIS VALUE HAS ALREADY BEEN COMPUTED.                  
            IF (IND .EQ. 6) GO TO 105                                   
               IF (IND .EQ. 5) THEN                                     
                  DO 102 K = 1,N                                        
                     W(K,1) = W(K,9)                                    
  102             CONTINUE                                              
               ELSE
                  CALL FCN(N,X,Y,W(1,1))
                  C(24) = C(24) + 1
               ENDIF                                                    
  105       CONTINUE                                                    
C                                                                       
C    CALCULATE HMIN - USE DEFAULT UNLESS VALUE PRESCRIBED               
       C(13) = C(3)                                                     
       IF (C(3) .NE. 0.D0) GO TO 165                                    
C         CALCULATE DEFAULT VALUE OF HMIN                               
C         FIRST CALCULATE WEIGHTED NORM Y - C(12) - AS SPECIFIED        
C         BY THE ERROR CONTROL INDICATOR C(1)                           
          TEMP = 0.D0                                                   
          IF (C(1) .NE. 1.D0) GO TO 115                                 
C  ABSOLUTE ERROR CONTROL - WEIGHTS ARE 1                               
          DO 110 K = 1, N                                               
             TEMP = DMAX1(TEMP, DABS(Y(K)))                             
  110     CONTINUE                                                      
          C(12) = TEMP                                                  
          GO TO 160                                                     
  115  IF (C(1) .NE. 2.D0) GO TO 120                                    
C  RELATIVE ERROR CONTROL - WEIGHTS ARE 1/DABS(Y(K)) SO                 
C  WEIGHTED NORM Y IS 1                                                 
          C(12) = 1.D0                                                  
          GO TO 160                                                     
  120  IF (C(1) .NE. 3.D0) GO TO 130                                    
C  WEIGHTS ARE 1/MAX(C(2),ABS(Y(K)))                                    
          DO 125 K = 1, N                                               
             TEMP = DMAX1(TEMP, DABS(Y(K))/C(2))                        
  125     CONTINUE                                                      
          C(12) = DMIN1(TEMP, 1.D0)                                     
          GO TO 160                                                     
  130  IF (C(1) .NE. 4.D0) GO TO 140                                    
C  WEIGHTS ARE 1/MAX(C(K+30),ABS(Y(K)))                                 
          DO 135 K = 1, N                                               
             TEMP = DMAX1(TEMP, DABS(Y(K))/C(K+30))                     
  135     CONTINUE                                                      
          C(12) = DMIN1(TEMP, 1.D0)                                     
          GO TO 160                                                     
  140  IF (C(1) .NE. 5.D0) GO TO 150                                    
C  WEIGHTS ARE 1/C(K+30)                                                
          DO 145 K = 1, N                                               
             TEMP = DMAX1(TEMP, DABS(Y(K))/C(K+30))                     
  145     CONTINUE                                                      
          C(12) = TEMP                                                  
          GO TO 160                                                     
  150  CONTINUE                                                         
C  DEFAULT CASE - WEIGHTS ARE 1/MAX(1,ABS(Y(K)))                        
          DO 155 K = 1, N                                               
             TEMP = DMAX1(TEMP, DABS(Y(K)))                             
  155     CONTINUE                                                      
          C(12) = DMIN1(TEMP, 1.D0)                                     
  160  CONTINUE                                                         
       C(13) = 10.D0*DMAX1(C(11),C(10)*DMAX1(C(12)/TOL,DABS(X)))        
  165  CONTINUE                                                         
C                                                                       
C    CALCULATE SCALE - USE DEFAULT UNLESS VALUE PRESCRIBED              
       C(15) = C(5)                                                     
       IF (C(5) .EQ. 0.D0) C(15) = 1.D0                                 
C                                                                       
C    CALCULATE HMAX - CONSIDER 4 CASES                                  
C    CASE 1 BOTH HMAX AND SCALE PRESCRIBED                              
       IF (C(6).NE.0.D0 .AND. C(5).NE.0.D0)                             
     +  C(16) = DMIN1(C(6), 2.D0/C(5))                                  
C    CASE 2 - HMAX PRESCRIBED, BUT SCALE NOT                            
       IF (C(6).NE.0.D0 .AND. C(5).EQ.0.D0) C(16) = C(6)                
C    CASE 3 - HMAX NOT PRESCRIBED, BUT SCALE IS                         
       IF (C(6).EQ.0.D0 .AND. C(5).NE.0.D0) C(16) = 2.D0/C(5)           
C    CASE 4 - NEITHER HMAX NOR SCALE IS PROVIDED                        
       IF (C(6).EQ.0.D0 .AND. C(5).EQ.0.D0) C(16) = 2.D0                
C                                                                       
C***********ERROR RETURN (WITH IND=-2) IF HMIN .GT. HMAX                
       IF (C(13) .LE. C(16)) GO TO 170                                  
          IND = -2                                                      
          RETURN                                                        
  170  CONTINUE                                                         
C                                                                       
C    CALCULATE PRELIMINARY HMAG - CONSIDER 3 CASES                      
       IF (IND .GT. 2) GO TO 175                                        
C    CASE 1 - INITIAL ENTRY - USE PRESCRIBED VALUE OF HSTART, IF        
C       ANY, ELSE DEFAULT                                               
          C(14) = C(4)                                                  
          IF (C(4) .EQ. 0.D0) C(14) = C(16)*TOL**(1./6.)                
          GO TO 185                                                     
  175  IF (C(23) .GT. 1.D0) GO TO 180                                   
C    CASE 2 - AFTER A SUCCESSFUL STEP, OR AT MOST  ONE  FAILURE,        
C       USE MIN(2, .9*(TOL/EST)**(1/6))*HMAG, BUT AVOID POSSIBLE        
C       OVERFLOW. THEN AVOID REDUCTION BY MORE THAN HALF.               
          TEMP = 2.D0*C(14)                                             
          IF (TOL .LT. (2.D0/.9D0)**6*C(19))                            
     +          TEMP = .9D0*(TOL/C(19))**(1./6.)*C(14)                  
          C(14) = DMAX1(TEMP, .5D0*C(14))                               
          GO TO 185                                                     
  180  CONTINUE                                                         
C    CASE 3 - AFTER TWO OR MORE SUCCESSIVE FAILURES                     
          C(14) = .5D0*C(14)                                            
  185  CONTINUE                                                         
C                                                                       
C    CHECK AGAINST HMAX                                                 
       C(14) = DMIN1(C(14), C(16))                                      
C                                                                       
C    CHECK AGAINST HMIN                                                 
       C(14) = DMAX1(C(14), C(13))                                      
C                                                                       
C***********INTERRUPT NO 1 (WITH IND=4) IF REQUESTED                    
       IF (C(8) .EQ. 0.D0) GO TO 1111                                   
          IND = 4                                                       
          RETURN                                                        
C    RESUME HERE ON RE-ENTRY WITH IND .EQ. 4   ........RE-ENTRY..       
 1111  CONTINUE                                                         
C                                                                       
C    CALCULATE HMAG, XTRIAL - DEPENDING ON PRELIMINARY HMAG, XEND       
       IF (C(14) .GE. DABS(XEND - X)) GO TO 190                         
C       DO NOT STEP MORE THAN HALF WAY TO XEND                          
          C(14) = DMIN1(C(14), .5D0*DABS(XEND - X))                     
          C(17) = X + DSIGN(C(14), XEND - X)                            
          GO TO 195                                                     
  190  CONTINUE                                                         
C       HIT XEND EXACTLY                                                
          C(14) = DABS(XEND - X)                                        
          C(17) = XEND                                                  
  195  CONTINUE                                                         
C                                                                       
C    CALCULATE HTRIAL                                                   
       C(18) = C(17) - X                                                
C                                                                       
C END STAGE 1                                                           
C                                                                       
C ***************************************************************       
C * STAGE 2 - CALCULATE YTRIAL (ADDING 10 TO NO OF  FCN  EVALS).*       
C * W(*,2), ... W(*,12) HOLD  INTERMEDIATE  RESULTS  NEEDED  IN *       
C * STAGE 3. W(*,13) IS TEMPORARY STORAGE UNTIL FINALLY IT HOLDS*       
C * YTRIAL.       *                                                     
C ***************************************************************       
C                                                                       
      TEMP = C(18)/1398169080000.D0                                     
C                                                                       
      DO 200 K = 1, N                                                   
         W(K,13) = Y(K) + TEMP*W(K,1)*233028180000.D0                   
  200 CONTINUE
      CALL FCN(N, X + C(18)/6.D0, W(1,13), W(1,2))                      
C
      DO 205 K = 1, N                                                   
         W(K,13) = Y(K) + TEMP*(W(K,1)*74569017600.D0                   
     +                       + W(K,2)*298276070400.D0)                  
  205 CONTINUE
      CALL FCN(N, X + C(18)*(4.D0/15.D0), W(1,13), W(1,3))              
C
      DO 210 K = 1, N                                                   
         W(K,13) = Y(K) + TEMP*(W(K,1)*1165140900000.D0                 
     +                       - W(K,2)*3728450880000.D0                  
     +                       + W(K,3)*3495422700000.D0 )                
  210 CONTINUE
      CALL FCN(N, X + C(18)*(2.D0/3.D0), W(1,13), W(1,4))               
C
      DO 215 K = 1, N                                                   
         W(K,13) = Y(K) + TEMP*( - W(K,1)*3604654659375.D0              
     +                          + W(K,2)*12816549900000.D0              
     +                          - W(K,3)*9284716546875.D0               
     +                          + W(K,4)*1237962206250.D0 )             
  215 CONTINUE
      CALL FCN(N, X + C(18)*(5.D0/6.D0), W(1,13), W(1,5))               
C                                                                       
C      DO 220 K = 1, N                                                   
C         W(K,13) = Y(K) + TEMP*(W(K,1)*3355605792000.D0                 
C     +                       - W(K,2)*11185352640000.D0                 
C     +                       + W(K,3)*9172628850000.D0                  
C     +                       - W(K,4)*427218330000.D0                   
C     +                       + W(K,5)*482505408000.D0)                  
C  220 CONTINUE
C      CALL FCN(N, X + C(18), W(1,13), W(1,6))                           
C                                                                       
      DO 225 K = 1, N                                                   
         W(K,13) = Y(K) + TEMP*( - W(K,1)*770204740536.D0               
     +                          + W(K,2)*2311639545600.D0               
     +                          - W(K,3)*1322092233000.D0               
     +                          - W(K,4)*453006781920.D0                
     +                          + W(K,5)*326875481856.D0)               
  225 CONTINUE
      CALL FCN(N, X + C(18)/15.D0, W(1,13), W(1,7))                     
C                                                                       
      DO 230 K = 1, N                                                   
         W(K,13) = Y(K) + TEMP*(W(K,1)*2845924389000.D0                 
     +                       - W(K,2)*9754668000000.D0                  
     +                       + W(K,3)*7897110375000.D0                  
     +                       - W(K,4)*192082660000.D0                   
     +                       + W(K,5)*400298976000.D0                   
     +                       + W(K,7)*201586000000.D0)                  
  230 CONTINUE
      CALL FCN(N, X + C(18), W(1,13), W(1,8))                           
C                                                                       
C    CALCULATE YTRIAL, THE EXTRAPOLATED APPROXIMATION AND STORE         
C       IN W(*,13)                                                      
      DO 235 K = 1, N                                                   
         W(K,13) = Y(K) + TEMP*(W(K,1)*104862681000.D0                  
     +                       + W(K,3)*545186250000.D0                   
     +                       + W(K,4)*446637345000.D0                   
     +                       + W(K,5)*188806464000.D0                   
     +                       + W(K,7)*15076875000.D0                    
     +                       + W(K,8)*97599465000.D0)                   
  235 CONTINUE                                                          
C
C NOW CALCULATE K9,K10 AND Y5 (STORED IN W(*,2))                        
C
      CALL FCN (N,  C(17), W(1,13), W(1,9) )                            
C                                                                       
      DO 240 I = 1,N                                                    
         W(I,2) = Y(I) + C(18) * ( A1*W(I,1) + A3*W(I,3) + A4*W(I,4)    
     +                  + A5*W(I,5) + A7*W(I,7) + A8*W(I,8) 
     +                  + A0*W(I,9) )                                   
 240  CONTINUE                                                          
C
      CALL FCN (N,  X + .5D0*C(18), W(1,2), W(1,10) )                   
C
CNOW USE THESE VALUES TO CALCULATE K11 AND K12 USING                    
C       W(*,14) AS TEMPORARY STROAGE.                                   
C                                                                       
      DO 250 I = 1,N                                                    
         W(I,14) = BY0*Y(I) + BY1*W(I,13) + BY5*W(I,2)                  
     +             + C(18) * ( BF0*W(I,1) + BF1*W(I,9) + BF5*W(I,10))   
 250  CONTINUE                                                          
C
      CALL FCN (N,  X + .90D0*C(18), W(1,14), W(1,11) )                 
C                                                                       
      DO 260 I = 1,N                                                    
         W(I,14) = CY0*Y(I) + CY1*W(I,13) + CY5*W(I,2)                  
     +          + C(18) * ( CF0*W(I,1) + CF1*W(I,9) + CF5*W(I,10))      
 260  CONTINUE                                                          
C
      CALL FCN (N,  X + .95D0*C(18), W(1,14), W(1,12) )                 
C
C    ADD 11 TO THE NO OF FCN EVALS                                      
      C(24) = C(24) + 10.D0                                             
C                                                                       
C END STAGE 2                                                           
C                                                                       
C        ***************************************************************
C        * STAGE 3 - CALCULATE THE ERROR ESTIMATE EST. FIRST CALCULATE *
C        * THE UNWEIGHTED  ABSOLUTE  ERROR  ESTIMATE VECTOR AND STORE  *
C        * IT IN W(*,2). THEN CALCULATE THE WEIGHTED MAX NORM OF W(*,2)*
C        * AS SPECIFIED BY THE ERROR CONTROL INDICATOR C(1).           *
C        ***************************************************************
C                                                                       
C           CALCULATE THE DEFECT AT X + .5*H AND USE IT AS AN ESTIMA
C           OF THE MAXIMUM DEFECT ON THIS STEP.                         
            TEMP = .1943D0                                             
            CALL INTRP(N,X,Y,X + TEMP*C(18),W(1,2),C(18),NW,W)
            CALL FCN(N,X + TEMP*C(18),W(1,2),W(1,14))
            CALL DERIV(N,X,Y,X + TEMP*C(18),W(1,2),C(18),NW,W)          
            DO 300 K = 1,N
               W(K,2) = W(K,2) - W(K,14)                                
  300       CONTINUE                                                    
            C(24) = C(24) + 1.D0                                        
C                                                                       
C    CALCULATE THE WEIGHTED MAX NORM OF W(*,2) AS SPECIFIED BY          
C    THE ERROR CONTROL INDICATOR C(1)                                   
      TEMP = 0.D0                                                       
      IF (C(1) .NE. 1.D0) GO TO 310                                     
C       ABSOLUTE ERROR CONTROL                                          
         DO 305 K = 1, N                                                
            TEMP = DMAX1(TEMP,DABS(W(K,2)))                             
  305    CONTINUE                                                       
         GO TO 360                                                      
  310 IF (C(1) .NE. 2.D0) GO TO 320                                     
C       RELATIVE ERROR CONTROL                                          
         DO 315 K = 1, N                                                
            TEMP = DMAX1(TEMP, DABS(W(K,2)/Y(K)))                       
  315    CONTINUE                                                       
         GO TO 360                                                      
  320 IF (C(1) .NE. 3.D0) GO TO 330                                     
C       WEIGHTS ARE 1/MAX(C(2),ABS(Y(K)))                               
         DO 325 K = 1, N                                                
            TEMP = DMAX1(TEMP, DABS(W(K,2))                             
     +                    / DMAX1(C(2), DABS(Y(K))) )                   
  325    CONTINUE                                                       
         GO TO 360                                                      
  330 IF (C(1) .NE. 4.D0) GO TO 340                                     
C       WEIGHTS ARE 1/MAX(C(K+30),ABS(Y(K)))                            
         DO 335 K = 1, N                                                
            TEMP = DMAX1(TEMP, DABS(W(K,2))                             
     +                    / DMAX1(C(K+30), DABS(Y(K))) )                
  335    CONTINUE                                                       
         GO TO 360                                                      
  340 IF (C(1) .NE. 5.D0) GO TO 350                                     
C       WEIGHTS ARE 1/C(K+30)                                           
         DO 345 K = 1, N                                                
            TEMP = DMAX1(TEMP, DABS(W(K,2)/C(K+30)))                    
  345    CONTINUE                                                       
         GO TO 360                                                      
  350 CONTINUE                                                          
C       DEFAULT CASE - WEIGHTS ARE 1/MAX(1,ABS(Y(K)))                   
         DO 355 K = 1, N                                                
            TEMP = DMAX1(TEMP, DABS(W(K,2))                             
     +                     / DMAX1(1.D0, DABS(Y(K))) )                  
  355    CONTINUE                                                       
  360 CONTINUE                                                          
C                                                                       
C    CALCULATE EST - (THE WEIGHTED MAX NORM OF W(*,2))*HMAG*SCALE       
      C(19) = TEMP                                                      
C                                                                       
C END STAGE 3                                                           
C                                                                       
C ***************************************************************       
C * STAGE 4 - MAKE DECISIONS.       *                                   
C ***************************************************************       
C                                                                       
C    SET IND=5 IF STEP ACCEPTABLE, ELSE SET IND=6                       
      IND = 5                                                           
      IF (C(19) .GT. TOL) IND = 6                                       
C                                                                       
C***********INTERRUPT NO 2 IF REQUESTED                                 
      IF (C(9) .EQ. 0.D0) GO TO 2222                                    
         RETURN                                                         
C    RESUME HERE ON RE-ENTRY WITH IND .EQ. 5 OR 6   ...RE-ENTRY..       
 2222 CONTINUE                                                          
C                                                                       
      IF (IND .EQ. 6) GO TO 410                                         
C       STEP ACCEPTED (IND .EQ. 5), SO UPDATE X, Y FROM XTRIAL,         
C  YTRIAL, ADD 1 TO THE NO OF SUCCESSFUL STEPS, AND SET                 
C  THE NO OF SUCCESSIVE FAILURES TO ZERO                                
         X = C(17)                                                      
         DO 400 K = 1, N                                                
            Y(K) = W(K,13)                                              
  400    CONTINUE                                                       
         C(22) = C(22) + 1.D0                                           
         C(23) = 0.D0                                                   
C**************RETURN(WITH IND=3, XEND SAVED, FLAG SET) IF X .EQ. XEND  
         IF (X .NE. XEND) GO TO 405                                     
            IND = 3                                                     
            C(20) = XEND                                                
            C(21) = 1.D0                                                
            RETURN                                                      
  405    CONTINUE                                                       
            GO TO 420                                                   
  410 CONTINUE                                                          
C       STEP NOT ACCEPTED (IND .EQ. 6), SO ADD 1 TO THE NO OF           
C  SUCCESSIVE FAILURES                                                  
         C(23) = C(23) + 1.D0                                           
C**************ERROR RETURN (WITH IND=-3) IF HMAG .LE. HMIN             
         IF (C(14) .GT. C(13)) GO TO 415                                
            IND = -3                                                    
            RETURN                                                      
  415    CONTINUE                                                       
  420 CONTINUE                                                          
C                                                                       
C END STAGE 4                                                           
C                                                                       
      GO TO 99999                                                       
C     END LOOP                                                          
C                                                                       
C  BEGIN ABORT ACTION                                                   
  500 CONTINUE                                                          
C                                                                       
      WRITE(6,505) IND, TOL, X, N, C(13), XEND, NW, C(16), C(20),       
     +    C(22), C(23), C(24), (Y(K), K = 1, N)                         
  505 FORMAT( /// 1H0, 58HCOMPUTATION STOPPED IN DVERK WITH THE FOLLOWIN
     +G VALUES -                                                        
     + / 1H0, 5HIND =, I4, 5X, 6HTOL  =, 1PD13.6, 5X, 11HX         =,   
     +1PD22.15                                                          
     + / 1H , 5HN   =, I4, 5X, 6HHMIN =, 1PD13.6, 5X, 11HXEND      =,   
     +1PD22.15                                                          
     + / 1H , 5HNW  =, I4, 5X, 6HHMAX =, 1PD13.6, 5X, 11HPREV XEND =,   
     +1PD22.15                                                          
     + / 1H0, 14X, 27HNO OF SUCCESSFUL STEPS  =  , 0PF8.0               
     + / 1H , 14X, 27HNO OF SUCCESSIVE FAILURES =, 0PF8.0               
     + / 1H , 14X, 27HNO OF FUNCTION EVALS  =    , 0PF8.0               
     + / 1H0, 23HTHE COMPONENTS OF Y ARE                                
     + // (1H , 1P5D24.15)       )                                      
C                                                                       
      STOP                                                              
C                                                                       
C  END ABORT ACTION                                                     
C                                                                       
      END                                                               
C                                                                       
