      SUBROUTINE DVERK(N, FCN, X, Y, XEND, TOL, IND, C, NW, W)
      INTEGER N, IND, NW, K
      DOUBLE PRECISION X, Y(N), XEND, TOL, C(*), W(NW,9), TEMP
C Added external statement for fcn to avoid a warning message.
      external fcn
C
C***********************************************************************
C								       *
C     PURPOSE - THIS IS A RUNGE-KUTTA  SUBROUTINE  BASED  ON  VERNER'S *
C FIFTH AND SIXTH ORDER PAIR OF FORMULAS FOR FINDING APPROXIMATIONS TO *
C THE SOLUTION OF  A  SYSTEM  OF  FIRST  ORDER	ORDINARY  DIFFERENTIAL *
C EQUATIONS  WITH  INITIAL  CONDITIONS. IT ATTEMPTS TO KEEP THE GLOBAL *
C ERROR PROPORTIONAL TO  A  TOLERANCE  SPECIFIED  BY  THE  USER.  (THE *
C PROPORTIONALITY  DEPENDS  ON THE KIND OF ERROR CONTROL THAT IS USED, *
C AS WELL AS THE DIFFERENTIAL EQUATION AND THE RANGE OF INTEGRATION.)  *
C								       *
C     VARIOUS OPTIONS ARE AVAILABLE TO THE USER,  INCLUDING  DIFFERENT *
C KINDS  OF  ERROR CONTROL, RESTRICTIONS ON STEP SIZES, AND INTERRUPTS *
C WHICH PERMIT THE USER TO EXAMINE THE STATE OF THE  CALCULATION  (AND *
C PERHAPS MAKE MODIFICATIONS) DURING INTERMEDIATE STAGES.	       *
C								       *
C     THE PROGRAM IS EFFICIENT FOR NON-STIFF SYSTEMS.  HOWEVER, A GOOD *
C VARIABLE-ORDER-ADAMS	METHOD	WILL PROBABLY BE MORE EFFICIENT IF THE *
C FUNCTION EVALUATIONS ARE VERY COSTLY.  SUCH A METHOD WOULD  ALSO  BE *
C MORE SUITABLE IF ONE WANTED TO OBTAIN A LARGE NUMBER OF INTERMEDIATE *
C SOLUTION VALUES BY INTERPOLATION, AS MIGHT BE THE CASE  FOR  EXAMPLE *
C WITH GRAPHICAL OUTPUT.					       *
C								       *
C				     HULL-ENRIGHT-JACKSON   1/10/76    *
C								       *
C***********************************************************************
C								       *
C     USE - THE USER MUST SPECIFY EACH OF THE FOLLOWING 	       *
C								       *
C     N  NUMBER OF EQUATIONS					       *
C								       *
C   FCN  NAME OF SUBROUTINE FOR EVALUATING FUNCTIONS - THE  SUBROUTINE *
C	    ITSELF MUST ALSO BE PROVIDED BY THE USER - IT SHOULD BE OF *
C	    THE FOLLOWING FORM					       *
C	       SUBROUTINE FCN(N, X, Y, YPRIME)			       *
C	       INTEGER N					       *
C	       DOUBLE PRECISION X, Y(N), YPRIME(N)		       *
C		       *** ETC ***				       *
C	    AND IT SHOULD EVALUATE YPRIME, GIVEN N, X AND Y	       *
C								       *
C     X  INDEPENDENT VARIABLE - INITIAL VALUE SUPPLIED BY USER	       *
C								       *
C     Y  DEPENDENT VARIABLE - INITIAL VALUES OF COMPONENTS Y(1), Y(2), *
C	    ..., Y(N) SUPPLIED BY USER				       *
C								       *
C  XEND  VALUE OF X TO WHICH INTEGRATION IS TO BE CARRIED OUT - IT MAY *
C	    BE LESS THAN THE INITIAL VALUE OF X 		       *
C								       *
C   TOL  TOLERANCE - THE SUBROUTINE ATTEMPTS TO CONTROL A NORM OF  THE *
C	    LOCAL  ERROR  IN  SUCH  A  WAY  THAT  THE  GLOBAL ERROR IS *
C	    PROPORTIONAL TO TOL. IN SOME PROBLEMS THERE WILL BE ENOUGH *
C	    DAMPING  OF  ERRORS, AS WELL AS SOME CANCELLATION, SO THAT *
C	    THE GLOBAL ERROR WILL BE LESS THAN TOL. ALTERNATIVELY, THE *
C	    CONTROL   CAN   BE	VIEWED	AS  ATTEMPTING	TO  PROVIDE  A *
C	    CALCULATED VALUE OF Y AT XEND WHICH IS THE EXACT  SOLUTION *
C	    TO	THE  PROBLEM Y' = F(X,Y) + E(X) WHERE THE NORM OF E(X) *
C	    IS PROPORTIONAL TO TOL.  (THE NORM	IS  A  MAX  NORM  WITH *
C	    WEIGHTS  THAT  DEPEND ON THE ERROR CONTROL STRATEGY CHOSEN *
C	    BY THE USER.  THE DEFAULT WEIGHT FOR THE K-TH COMPONENT IS *
C	    1/MAX(1,ABS(Y(K))),  WHICH THEREFORE PROVIDES A MIXTURE OF *
C	    ABSOLUTE AND RELATIVE ERROR CONTROL.)		       *
C								       *
C   IND  INDICATOR - ON INITIAL ENTRY IND MUST BE SET EQUAL TO	EITHER *
C	    1  OR  2. IF THE USER DOES NOT WISH TO USE ANY OPTIONS, HE *
C	    SHOULD SET IND TO 1 - ALL THAT REMAINS FOR THE USER TO  DO *
C	    THEN  IS  TO  DECLARE C AND W, AND TO SPECIFY NW. THE USER *
C	    MAY ALSO  SELECT  VARIOUS  OPTIONS	ON  INITIAL  ENTRY  BY *
C	    SETTING IND = 2 AND INITIALIZING THE FIRST 9 COMPONENTS OF *
C	    C AS DESCRIBED IN THE NEXT SECTION.  HE MAY ALSO  RE-ENTER *
C	    THE  SUBROUTINE  WITH IND = 3 AS MENTIONED AGAIN BELOW. IN *
C	    ANY EVENT, THE SUBROUTINE RETURNS WITH IND EQUAL TO        *
C	       3 AFTER A NORMAL RETURN				       *
C	       4, 5, OR 6 AFTER AN INTERRUPT (SEE OPTIONS C(8), C(9))  *
C	       -1, -2, OR -3 AFTER AN ERROR CONDITION (SEE BELOW)      *
C								       *
C     C  COMMUNICATIONS VECTOR - THE DIMENSION MUST BE GREATER THAN OR *
C	    EQUAL TO 24, UNLESS OPTION C(1) = 4 OR 5 IS USED, IN WHICH *
C	    CASE THE DIMENSION MUST BE GREATER THAN OR EQUAL TO N+30   *
C								       *
C    NW  FIRST DIMENSION OF WORKSPACE W -  MUST  BE  GREATER  THAN  OR *
C	    EQUAL TO N						       *
C								       *
C     W  WORKSPACE MATRIX - FIRST DIMENSION MUST BE NW AND SECOND MUST *
C	    BE GREATER THAN OR EQUAL TO 9			       *
C								       *
C     THE SUBROUTINE  WILL  NORMALLY  RETURN  WITH  IND  =  3,	HAVING *
C REPLACED THE INITIAL VALUES OF X AND Y WITH, RESPECTIVELY, THE VALUE *
C OF XEND AND AN APPROXIMATION TO Y AT XEND.  THE  SUBROUTINE  CAN  BE *
C CALLED  REPEATEDLY  WITH NEW VALUES OF XEND WITHOUT HAVING TO CHANGE *
C ANY OTHER ARGUMENT.  HOWEVER, CHANGES IN TOL, OR ANY OF THE  OPTIONS *
C DESCRIBED BELOW, MAY ALSO BE MADE ON SUCH A RE-ENTRY IF DESIRED.     *
C								       *
C     THREE ERROR RETURNS ARE ALSO POSSIBLE, IN WHICH  CASE  X	AND  Y *
C WILL BE THE MOST RECENTLY ACCEPTED VALUES -			       *
C     WITH IND = -3 THE SUBROUTINE WAS UNABLE  TO  SATISFY  THE  ERROR *
C	 REQUIREMENT  WITH A PARTICULAR STEP-SIZE THAT IS LESS THAN OR *
C	 EQUAL TO HMIN, WHICH MAY MEAN THAT TOL IS TOO SMALL	       *
C     WITH IND = -2 THE VALUE OF HMIN  IS  GREATER  THAN  HMAX,  WHICH *
C	 PROBABLY  MEANS  THAT THE REQUESTED TOL (WHICH IS USED IN THE *
C	 CALCULATION OF HMIN) IS TOO SMALL			       *
C     WITH IND = -1 THE ALLOWED MAXIMUM NUMBER OF FCN EVALUATIONS  HAS *
C	 BEEN  EXCEEDED,  BUT  THIS  CAN ONLY OCCUR IF OPTION C(7), AS *
C	 DESCRIBED IN THE NEXT SECTION, HAS BEEN USED		       *
C								       *
C     THERE ARE SEVERAL CIRCUMSTANCES THAT WILL CAUSE THE CALCULATIONS *
C TO  BE  TERMINATED,  ALONG WITH OUTPUT OF INFORMATION THAT WILL HELP *
C THE USER DETERMINE THE CAUSE OF  THE	TROUBLE.  THESE  CIRCUMSTANCES *
C INVOLVE  ENTRY WITH ILLEGAL OR INCONSISTENT VALUES OF THE ARGUMENTS, *
C SUCH AS ATTEMPTING A NORMAL  RE-ENTRY  WITHOUT  FIRST  CHANGING  THE *
C VALUE OF XEND, OR ATTEMPTING TO RE-ENTER WITH IND LESS THAN ZERO.    *
C								       *
C***********************************************************************
C								       *
C     OPTIONS - IF THE SUBROUTINE IS ENTERED WITH IND = 1, THE FIRST 9 *
C COMPONENTS OF THE COMMUNICATIONS VECTOR ARE INITIALIZED TO ZERO, AND *
C THE SUBROUTINE USES ONLY DEFAULT VALUES  FOR	EACH  OPTION.  IF  THE *
C SUBROUTINE  IS  ENTERED  WITH IND = 2, THE USER MUST SPECIFY EACH OF *
C THESE 9 COMPONENTS - NORMALLY HE WOULD FIRST SET THEM ALL  TO  ZERO, *
C AND  THEN  MAKE  NON-ZERO  THOSE  THAT  CORRESPOND TO THE PARTICULAR *
C OPTIONS HE WISHES TO SELECT. IN ANY EVENT, OPTIONS MAY BE CHANGED ON *
C RE-ENTRY  TO	THE  SUBROUTINE  -  BUT IF THE USER CHANGES ANY OF THE *
C OPTIONS, OR TOL, IN THE COURSE OF A CALCULATION HE SHOULD BE CAREFUL *
C ABOUT  HOW  SUCH CHANGES AFFECT THE SUBROUTINE - IT MAY BE BETTER TO *
C RESTART WITH IND = 1 OR 2. (COMPONENTS 10 TO 24 OF C ARE USED BY THE *
C PROGRAM  -  THE INFORMATION IS AVAILABLE TO THE USER, BUT SHOULD NOT *
C NORMALLY BE CHANGED BY HIM.)					       *
C								       *
C  C(1)  ERROR CONTROL INDICATOR - THE NORM OF THE LOCAL ERROR IS  THE *
C	    MAX  NORM  OF  THE	WEIGHTED  ERROR  ESTIMATE  VECTOR, THE *
C	    WEIGHTS BEING DETERMINED ACCORDING TO THE VALUE OF C(1) -  *
C	       IF C(1)=1 THE WEIGHTS ARE 1 (ABSOLUTE ERROR CONTROL)    *
C	       IF C(1)=2 THE WEIGHTS ARE 1/ABS(Y(K))  (RELATIVE  ERROR *
C		  CONTROL)					       *
C	       IF C(1)=3 THE  WEIGHTS  ARE  1/MAX(ABS(C(2)),ABS(Y(K))) *
C		  (RELATIVE  ERROR  CONTROL,  UNLESS ABS(Y(K)) IS LESS *
C		  THAN THE FLOOR VALUE, ABS(C(2)) )		       *
C	       IF C(1)=4 THE WEIGHTS ARE 1/MAX(ABS(C(K+30)),ABS(Y(K))) *
C		  (HERE INDIVIDUAL FLOOR VALUES ARE USED)	       *
C	       IF C(1)=5 THE WEIGHTS ARE 1/ABS(C(K+30)) 	       *
C	       FOR ALL OTHER VALUES OF C(1), INCLUDING	C(1) = 0,  THE *
C		  DEFAULT  VALUES  OF  THE  WEIGHTS  ARE  TAKEN  TO BE *
C		  1/MAX(1,ABS(Y(K))), AS MENTIONED EARLIER	       *
C	    (IN THE TWO CASES C(1) = 4 OR 5 THE USER MUST DECLARE  THE *
C	    DIMENSION OF C TO BE AT LEAST N+30 AND MUST INITIALIZE THE *
C	    COMPONENTS C(31), C(32), ..., C(N+30).)		       *
C								       *
C  C(2)  FLOOR VALUE - USED WHEN THE INDICATOR C(1) HAS THE VALUE 3    *
C								       *
C  C(3)  HMIN SPECIFICATION - IF NOT ZERO, THE SUBROUTINE CHOOSES HMIN *
C	    TO BE ABS(C(3)) - OTHERWISE IT USES THE DEFAULT VALUE      *
C	       10*MAX(DWARF,RREB*MAX(WEIGHTED NORM Y/TOL,ABS(X))),     *
C	    WHERE DWARF IS A VERY SMALL POSITIVE  MACHINE  NUMBER  AND *
C	    RREB IS THE RELATIVE ROUNDOFF ERROR BOUND		       *
C								       *
C  C(4)  HSTART SPECIFICATION - IF NOT ZERO, THE SUBROUTINE  WILL  USE *
C	    AN	INITIAL  HMAG EQUAL TO ABS(C(4)), EXCEPT OF COURSE FOR *
C	    THE RESTRICTIONS IMPOSED BY HMIN AND HMAX  -  OTHERWISE IT *
C	    USES THE DEFAULT VALUE OF HMAX*(TOL)**(1/6) 	       *
C								       *
C  C(5)  SCALE SPECIFICATION - THIS IS INTENDED TO BE A MEASURE OF THE *
C	    SCALE OF THE PROBLEM - LARGER VALUES OF SCALE TEND TO MAKE *
C	    THE METHOD MORE RELIABLE, FIRST  BY  POSSIBLY  RESTRICTING *
C	    HMAX  (AS  DESCRIBED  BELOW) AND SECOND, BY TIGHTENING THE *
C	    ACCEPTANCE REQUIREMENT - IF C(5) IS ZERO, A DEFAULT  VALUE *
C	    OF	1  IS  USED.  FOR  LINEAR  HOMOGENEOUS	PROBLEMS  WITH *
C	    CONSTANT COEFFICIENTS, AN APPROPRIATE VALUE FOR SCALE IS A *
C	    NORM  OF  THE  ASSOCIATED  MATRIX.	FOR OTHER PROBLEMS, AN *
C	    APPROXIMATION TO  AN  AVERAGE  VALUE  OF  A  NORM  OF  THE *
C	    JACOBIAN ALONG THE TRAJECTORY MAY BE APPROPRIATE	       *
C								       *
C  C(6)  HMAX SPECIFICATION - FOUR CASES ARE POSSIBLE		       *
C	    IF C(6).NE.0 AND C(5).NE.0, HMAX IS TAKEN TO BE	       *
C	       MIN(ABS(C(6)),2/ABS(C(5)))			       *
C	    IF C(6).NE.0 AND C(5).EQ.0, HMAX IS TAKEN TO BE  ABS(C(6)) *
C	    IF C(6).EQ.0 AND C(5).NE.0, HMAX IS TAKEN TO BE	       *
C	       2/ABS(C(5))					       *
C	    IF C(6).EQ.0 AND C(5).EQ.0, HMAX IS GIVEN A DEFAULT  VALUE *
C	       OF 2						       *
C								       *
C  C(7)  MAXIMUM NUMBER OF FUNCTION EVALUATIONS  -  IF	NOT  ZERO,  AN *
C	    ERROR  RETURN WITH IND = -1 WILL BE CAUSED WHEN THE NUMBER *
C	    OF FUNCTION EVALUATIONS EXCEEDS ABS(C(7))		       *
C								       *
C  C(8)  INTERRUPT NUMBER  1  -  IF  NOT  ZERO,  THE  SUBROUTINE  WILL *
C	    INTERRUPT	THE  CALCULATIONS  AFTER  IT  HAS  CHOSEN  ITS *
C	    PRELIMINARY VALUE OF HMAG, AND JUST BEFORE CHOOSING HTRIAL *
C	    AND  XTRIAL  IN  PREPARATION FOR TAKING A STEP (HTRIAL MAY *
C	    DIFFER FROM HMAG IN SIGN, AND MAY  REQUIRE	ADJUSTMENT  IF *
C	    XEND  IS  NEAR) - THE SUBROUTINE RETURNS WITH IND = 4, AND *
C	    WILL RESUME CALCULATION AT THE POINT  OF  INTERRUPTION  IF *
C	    RE-ENTERED WITH IND = 4				       *
C								       *
C  C(9)  INTERRUPT NUMBER  2  -  IF  NOT  ZERO,  THE  SUBROUTINE  WILL *
C	    INTERRUPT	THE  CALCULATIONS  IMMEDIATELY	AFTER  IT  HAS *
C	    DECIDED WHETHER OR NOT TO ACCEPT THE RESULT  OF  THE  MOST *
C	    RECENT  TRIAL STEP, WITH IND = 5 IF IT PLANS TO ACCEPT, OR *
C	    IND = 6 IF IT PLANS TO REJECT -  Y(*)  IS  THE  PREVIOUSLY *
C	    ACCEPTED  RESULT, WHILE W(*,9) IS THE NEWLY COMPUTED TRIAL *
C	    VALUE, AND W(*,2) IS THE UNWEIGHTED ERROR ESTIMATE VECTOR. *
C	    THE  SUBROUTINE  WILL  RESUME CALCULATIONS AT THE POINT OF *
C	    INTERRUPTION ON RE-ENTRY WITH IND = 5 OR 6. (THE USER  MAY *
C	    CHANGE IND IN THIS CASE IF HE WISHES, FOR EXAMPLE TO FORCE *
C	    ACCEPTANCE OF A STEP THAT WOULD OTHERWISE BE REJECTED,  OR *
C	    VICE VERSA. HE CAN ALSO RESTART WITH IND = 1 OR 2.)        *
C								       *
C***********************************************************************
C								       *
C  SUMMARY OF THE COMPONENTS OF THE COMMUNICATIONS VECTOR	       *
C								       *
C     PRESCRIBED AT THE OPTION	     DETERMINED BY THE PROGRAM	       *
C	    OF THE USER 					       *
C								       *
C				     C(10) RREB(REL ROUNDOFF ERR BND)  *
C     C(1) ERROR CONTROL INDICATOR   C(11) DWARF (VERY SMALL MACH NO)  *
C     C(2) FLOOR VALUE		     C(12) WEIGHTED NORM Y	       *
C     C(3) HMIN SPECIFICATION	     C(13) HMIN 		       *
C     C(4) HSTART SPECIFICATION      C(14) HMAG 		       *
C     C(5) SCALE SPECIFICATION	     C(15) SCALE		       *
C     C(6) HMAX SPECIFICATION	     C(16) HMAX 		       *
C     C(7) MAX NO OF FCN EVALS	     C(17) XTRIAL		       *
C     C(8) INTERRUPT NO 1	     C(18) HTRIAL		       *
C     C(9) INTERRUPT NO 2	     C(19) EST			       *
C				     C(20) PREVIOUS XEND	       *
C				     C(21) FLAG FOR XEND	       *
C				     C(22) NO OF SUCCESSFUL STEPS      *
C				     C(23) NO OF SUCCESSIVE FAILURES   *
C				     C(24) NO OF FCN EVALS	       *
C								       *
C  IF C(1) = 4 OR 5, C(31), C(32), ... C(N+30) ARE FLOOR VALUES        *
C								       *
C  RREB and DWARF are machine dependent constants currently set so     *
C  that they should be appropriate for most machines.  However, it may *
C  be appropriate to change them when this program is installed on a   *
C  new machine.                            K.R.J.   3 Oct 1991.        *
C								       *
C***********************************************************************
C								       *
C  AN OVERVIEW OF THE PROGRAM					       *
C								       *
C     BEGIN INITIALIZATION, PARAMETER CHECKING, INTERRUPT RE-ENTRIES   *
C  ......ABORT IF IND OUT OF RANGE 1 TO 6			       *
C  .	 CASES - INITIAL ENTRY, NORMAL RE-ENTRY, INTERRUPT RE-ENTRIES  *
C  .	 CASE 1 - INITIAL ENTRY (IND .EQ. 1 OR 2)		       *
C  V........ABORT IF N.GT.NW OR TOL.LE.0			       *
C  .	    IF INITIAL ENTRY WITHOUT OPTIONS (IND .EQ. 1)	       *
C  .	       SET C(1) TO C(9) EQUAL TO ZERO			       *
C  .	    ELSE INITIAL ENTRY WITH OPTIONS (IND .EQ. 2)	       *
C  .	       MAKE C(1) TO C(9) NON-NEGATIVE			       *
C  .	       MAKE FLOOR VALUES NON-NEGATIVE IF THEY ARE TO BE USED   *
C  .	    END IF						       *
C  .	    INITIALIZE RREB, DWARF, PREV XEND, FLAG, COUNTS	       *
C  .	 CASE 2 - NORMAL RE-ENTRY (IND .EQ. 3)			       *
C  .........ABORT IF XEND REACHED, AND EITHER X CHANGED OR XEND NOT    *
C  .	    RE-INITIALIZE FLAG					       *
C  .	 CASE 3 - RE-ENTRY FOLLOWING AN INTERRUPT (IND .EQ. 4 TO 6)    *
C  V	    TRANSFER CONTROL TO THE APPROPRIATE RE-ENTRY POINT.......  *
C  .	 END CASES						    .  *
C  .  END INITIALIZATION, ETC.					    .  *
C  .								    V  *
C  .  LOOP THROUGH THE FOLLOWING 4 STAGES, ONCE FOR EACH TRIAL STEP .  *
C  .	 STAGE 1 - PREPARE					    .  *
C***********ERROR RETURN (WITH IND=-1) IF NO OF FCN EVALS TOO GREAT .  *
C  .	    CALC SLOPE (ADDING 1 TO NO OF FCN EVALS) IF IND .NE. 6  .  *
C  .	    CALC HMIN, SCALE, HMAX				    .  *
C***********ERROR RETURN (WITH IND=-2) IF HMIN .GT. HMAX	    .  *
C  .	    CALC PRELIMINARY HMAG				    .  *
C***********INTERRUPT NO 1 (WITH IND=4) IF REQUESTED.......RE-ENTRY.V  *
C  .	    CALC HMAG, XTRIAL AND HTRIAL			    .  *
C  .	 END STAGE 1						    .  *
C  V	 STAGE 2 - CALC YTRIAL (ADDING 7 TO NO OF FCN EVALS)	    .  *
C  .	 STAGE 3 - CALC THE ERROR ESTIMATE			    .  *
C  .	 STAGE 4 - MAKE DECISIONS				    .  *
C  .	    SET IND=5 IF STEP ACCEPTABLE, ELSE SET IND=6	    .  *
C***********INTERRUPT NO 2 IF REQUESTED....................RE-ENTRY.V  *
C  .	    IF STEP ACCEPTED (IND .EQ. 5)			       *
C  .	       UPDATE X, Y FROM XTRIAL, YTRIAL			       *
C  .	       ADD 1 TO NO OF SUCCESSFUL STEPS			       *
C  .	       SET NO OF SUCCESSIVE FAILURES TO ZERO		       *
C**************RETURN(WITH IND=3, XEND SAVED, FLAG SET) IF X .EQ. XEND *
C  .	    ELSE STEP NOT ACCEPTED (IND .EQ. 6) 		       *
C  .	       ADD 1 TO NO OF SUCCESSIVE FAILURES		       *
C**************ERROR RETURN (WITH IND=-3) IF HMAG .LE. HMIN	       *
C  .	    END IF						       *
C  .	 END STAGE 4						       *
C  .  END LOOP							       *
C  .								       *
C  BEGIN ABORT ACTION						       *
C     OUTPUT APPROPRIATE  MESSAGE  ABOUT  STOPPING  THE  CALCULATIONS, *
C	 ALONG WITH VALUES OF IND, N, NW, TOL, HMIN,  HMAX,  X,  XEND, *
C	 PREVIOUS XEND,  NO OF	SUCCESSFUL  STEPS,  NO	OF  SUCCESSIVE *
C	 FAILURES, NO OF FCN EVALS, AND THE COMPONENTS OF Y	       *
C     STOP							       *
C  END ABORT ACTION						       *
C								       *
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  UP	WITH  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	    CALCULATE SLOPE (ADDING 1 TO NO OF FCN EVALS) IF IND .NE. 6
	    IF (IND .EQ. 6) GO TO 105
	       CALL FCN(N, X, Y, W(1,1))
	       C(24) = C(24) + 1.D0
  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 7 TO NO OF  FCN  EVALS). *
C	 * W(*,2), ... W(*,8)  HOLD  INTERMEDIATE  RESULTS  NEEDED  IN *
C	 * STAGE 3. W(*,9) IS TEMPORARY STORAGE UNTIL FINALLY IT HOLDS *
C	 * YTRIAL.						       *
C	 ***************************************************************
C
	    TEMP = C(18)/1398169080000.D0
C
	    DO 200 K = 1, N
	       W(K,9) = Y(K) + TEMP*W(K,1)*233028180000.D0
  200	    CONTINUE
	    CALL FCN(N, X + C(18)/6.D0, W(1,9), W(1,2))
C
	    DO 205 K = 1, N
	       W(K,9) = 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,9), W(1,3))
C
	    DO 210 K = 1, N
	       W(K,9) = 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,9), W(1,4))
C
	    DO 215 K = 1, N
	       W(K,9) = 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,9), W(1,5))
C
	    DO 220 K = 1, N
	       W(K,9) = Y(K) + TEMP*(	W(K,1)*3355605792000.D0
     +				      - W(K,2)*11185352640000.D0
     +				      + W(K,3)*9172628850000.D0
     +				      - W(K,4)*427218330000.D0
     +				      + W(K,5)*482505408000.D0	)
  220	    CONTINUE
	    CALL FCN(N, X + C(18), W(1,9), W(1,6))
C
	    DO 225 K = 1, N
	       W(K,9) = 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,9), W(1,7))
C
	    DO 230 K = 1, N
	       W(K,9) = 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,9), W(1,8))
C
C	    CALCULATE YTRIAL, THE EXTRAPOLATED APPROXIMATION AND STORE
C	       IN W(*,9)
	    DO 235 K = 1, N
	       W(K,9) = 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	    ADD 7 TO THE NO OF FCN EVALS
	    C(24) = C(24) + 7.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 (PER UNIT *
C	 * STEP) FOR THE UNEXTRAPOLATED APPROXIMATION AND STORE IT  IN *
C	 * W(*,2).  THEN  CALCULATE THE WEIGHTED MAX NORM OF W(*,2) AS *
C	 * SPECIFIED BY THE ERROR  CONTROL  INDICATOR  C(1).  FINALLY, *
C	 * MODIFY  THIS RESULT TO PRODUCE EST, THE ERROR ESTIMATE (PER *
C	 * UNIT STEP) FOR THE EXTRAPOLATED APPROXIMATION YTRIAL.       *
C	 ***************************************************************
C
C	    CALCULATE THE UNWEIGHTED ABSOLUTE ERROR ESTIMATE VECTOR
	    DO 300 K = 1, N
	       W(K,2) = (   W(K,1)*8738556750.D0
     +			  + W(K,3)*9735468750.D0
     +			  - W(K,4)*9709507500.D0
     +			  + W(K,5)*8582112000.D0
     +			  + W(K,6)*95329710000.D0
     +			  - W(K,7)*15076875000.D0
     +			  - W(K,8)*97599465000.D0)/1398169080000.D0
  300	    CONTINUE
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	       - EST IS INTENDED TO BE A MEASURE OF THE ERROR  PER  UNIT
C	       STEP IN YTRIAL
	    C(19) = TEMP*C(14)*C(15)
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,9)
  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
