SUBROUTINE DRKGS

PURPOSE
   TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL
   EQUATIONS WITH GIVEN INITIAL VALUES.

USAGE
   CALL DRKGS (PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)
   PARAMETERS FCT AND OUTP REQUIRE AN EXTERNAL STATEMENT.

DESCRIPTION OF PARAMETERS
   PRMT   - DOUBLE PRECISION INPUT AND OUTPUT VECTOR WITH
	    DIMENSION GREATER THAN OR EQUAL TO 5, WHICH
	    SPECIFIES THE PARAMETERS OF THE INTERVAL AND OF
	    ACCURACY AND WHICH SERVES FOR COMMUNICATION BETWEEN
	    OUTPUT SUBROUTINE (FURNISHED BY THE USER) AND
	    SUBROUTINE DRKGS. EXCEPT PRMT(5) THE COMPONENTS
	    ARE NOT DESTROYED BY SUBROUTINE DRKGS AND THEY ARE
   PRMT(1)- LOWER BOUND OF THE INTERVAL (INPUT),
   PRMT(2)- UPPER BOUND OF THE INTERVAL (INPUT),
   PRMT(3)- INITIAL INCREMENT OF THE INDEPENDENT VARIABLE
	    (INPUT),
   PRMT(4)- UPPER ERROR BOUND (INPUT). IF ABSOLUTE ERROR IS
	    GREATER THAN PRMT(4), INCREMENT GETS HALVED.
	    IF INCREMENT IS LESS THAN PRMT(3) AND ABSOLUTE
	    ERROR LESS THAN PRMT(4)/50, INCREMENT GETS DOUBLED.
	    THE USER MAY CHANGE PRMT(4) BY MEANS OF HIS
	    OUTPUT SUBROUTINE.
   PRMT(5)- NO INPUT PARAMETER. SUBROUTINE DRKGS INITIALIZES
	    PRMT(5)=0. IF THE USER WANTS TO TERMINATE
	    SUBROUTINE DRKGS AT ANY OUTPUT POINT, HE HAS TO
	    CHANGE PRMT(5) TO NON-ZERO BY MEANS OF SUBROUTINE
	    OUTP. FURTHER COMPONENTS OF VECTOR PRMT ARE
	    FEASIBLE IF ITS DIMENSION IS DEFINED GREATER
	    THAN 5. HOWEVER SUBROUTINE DRKGS DOES NOT REQUIRE
	    AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL
	    FOR HANDING RESULT VALUES TO THE MAIN PROGRAM
	    (CALLING DRKGS) WHICH ARE OBTAINED BY SPECIAL
	    MANIPULATIONS WITH OUTPUT DATA IN SUBROUTINE OUTP.
   Y	  - DOUBLE PRECISION INPUT VECTOR OF INITIAL VALUES
	    (DESTROYED). LATERON Y IS THE RESULTING VECTOR OF
	    DEPENDENT VARIABLES COMPUTED AT INTERMEDIATE
	    POINTS X.
   DERY   - DOUBLE PRECISION INPUT VECTOR OF ERROR WEIGHTS
	    (DESTROYED). THE SUM OF ITS COMPONENTS MUST BE
	    EQUAL TO 1. LATERON DERY IS THE VECTOR OF
	    DERIVATIVES, WHICH BELONG TO FUNCTION VALUES Y AT
	    INTERMEDIATE POINTS X.
   NDIM   - AN INPUT VALUE, WHICH SPECIFIES THE NUMBER OF
	    EQUATIONS IN THE SYSTEM.
   IHLF   - AN OUTPUT VALUE, WHICH SPECIFIES THE NUMBER OF
	    BISECTIONS OF THE INITIAL INCREMENT. IF IHLF GETS
	    GREATER THAN 10, SUBROUTINE DRKGS RETURNS WITH
	    ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM. ERROR
	    MESSAGE IHLF=12 OR IHLF=13 APPEARS IN CASE
	    PRMT(3)=0 OR IN CASE SIGN(PRMT(3)).NE.SIGN(PRMT(2)-
	    PRMT(1)) RESPECTIVELY.
   FCT	  - THE NAME OF AN EXTERNAL SUBROUTINE USED. THIS
	    SUBROUTINE COMPUTES THE RIGHT HAND SIDES DERY OF
	    THE SYSTEM TO GIVEN VALUES X AND Y. ITS PARAMETER
	    LIST MUST BE X,Y,DERY. SUBROUTINE FCT SHOULD
	    NOT DESTROY X AND Y.
   OUTP   - THE NAME OF AN EXTERNAL OUTPUT SUBROUTINE USED.
	    ITS PARAMETER LIST MUST BE X,Y,DERY,IHLF,NDIM,PRMT.
	    NONE OF THESE PARAMETERS (EXCEPT, IF NECESSARY,
	    PRMT(4),PRMT(5),...) SHOULD BE CHANGED BY
	    SUBROUTINE OUTP. IF PRMT(5) IS CHANGED TO NON-ZERO,
	    SUBROUTINE DRKGS IS TERMINATED.
   AUX	  - DOUBLE PRECISION AUXILIARY STORAGE ARRAY WITH 8
	    ROWS AND NDIM COLUMNS.

REMARKS
   THE PROCEDURE TERMINATES AND RETURNS TO CALLING PROGRAM, IF
   (1) MORE THAN 10 BISECTIONS OF THE INITIAL INCREMENT ARE
       NECESSARY TO GET SATISFACTORY ACCURACY (ERROR MESSAGE
       IHLF=11),
   (2) INITIAL INCREMENT IS EQUAL TO 0 OR HAS WRONG SIGN
       (ERROR MESSAGES IHLF=12 OR IHLF=13),
   (3) THE WHOLE INTEGRATION INTERVAL IS WORKED THROUGH,
   (4) SUBROUTINE OUTP HAS CHANGED PRMT(5) TO NON-ZERO.

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
   THE EXTERNAL SUBROUTINES FCT(X,Y,DERY) AND
   OUTP(X,Y,DERY,IHLF,NDIM,PRMT) MUST BE FURNISHED BY THE USER.

METHOD
   EVALUATION IS DONE BY MEANS OF FOURTH ORDER RUNGE-KUTTA
   FORMULAE IN THE MODIFICATION DUE TO GILL. ACCURACY IS
   TESTED COMPARING THE RESULTS OF THE PROCEDURE WITH SINGLE
   AND DOUBLE INCREMENT.
   SUBROUTINE DRKGS AUTOMATICALLY ADJUSTS THE INCREMENT DURING
   THE WHOLE COMPUTATION BY HALVING OR DOUBLING. IF MORE THAN
   10 BISECTIONS OF THE INCREMENT ARE NECESSARY TO GET
   SATISFACTORY ACCURACY, THE SUBROUTINE RETURNS WITH
   ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM.
   TO GET FULL FLEXIBILITY IN OUTPUT, AN OUTPUT SUBROUTINE
   MUST BE FURNISHED BY THE USER.
   FOR REFERENCE, SEE
   RALSTON/WILF, MATHEMATICAL METHODS FOR DIGITAL COMPUTERS,
   WILEY, NEW YORK/LONDON, 1960, PP.110-120.