C 
C SUBROUTINE 'HPCG' FROM THE IBM SCIENTIFIC SUBROUTINE PACKAGE.
C
C NOTE:  TO CONFORM WITH THE FORTRAN77 STANDARD, DUMMY ARRAY DIMENSIONS
C        (1) HAVE BEEN CHANGED TO (*).
C
C     ..................................................................
C
C        SUBROUTINE HPCG
C
C        PURPOSE
C           TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY GENERAL
C           DIFFERENTIAL EQUATIONS WITH GIVEN INITIAL VALUES.
C
C        USAGE
C           CALL HPCG (PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)
C           PARAMETERS FCT AND OUTP REQUIRE AN EXTERNAL STATEMENT.
C
C        DESCRIPTION OF PARAMETERS
C           PRMT   - AN INPUT AND OUTPUT VECTOR WITH DIMENSION GREATER
C                    OR EQUAL TO 5, WHICH SPECIFIES THE PARAMETERS OF
C                    THE INTERVAL AND OF ACCURACY AND WHICH SERVES FOR
C                    COMMUNICATION BETWEEN OUTPUT SUBROUTINE (FURNISHED
C                    BY THE USER) AND SUBROUTINE HPCG. EXCEPT PRMT(5)
C                    THE COMPONENTS ARE NOT DESTROYED BY SUBROUTINE
C                    HPCG AND THEY ARE
C           PRMT(1)- LOWER BOUND OF THE INTERVAL (INPUT),
C           PRMT(2)- UPPER BOUND OF THE INTERVAL (INPUT),
C           PRMT(3)- INITIAL INCREMENT OF THE INDEPENDENT VARIABLE
C                    (INPUT),
C           PRMT(4)- UPPER ERROR BOUND (INPUT). IF ABSOLUTE ERROR IS
C                    GREATER THAN PRMT(4), INCREMENT GETS HALVED.
C                    IF INCREMENT IS LESS THAN PRMT(3) AND ABSOLUTE
C                    ERROR LESS THAN PRMT(4)/50, INCREMENT GETS DOUBLED.
C                    THE USER MAY CHANGE PRMT(4) BY MEANS OF HIS
C                    OUTPUT SUBROUTINE.
C           PRMT(5)- NO INPUT PARAMETER. SUBROUTINE HPCG INITIALIZES
C                    PRMT(5)=0. IF THE USER WANTS TO TERMINATE
C                    SUBROUTINE HPCG AT ANY OUTPUT POINT, HE HAS TO
C                    CHANGE PRMT(5) TO NON-ZERO BY MEANS OF SUBROUTINE
C                    OUTP. FURTHER COMPONENTS OF VECTOR PRMT ARE
C                    FEASIBLE IF ITS DIMENSION IS DEFINED GREATER
C                    THAN 5. HOWEVER SUBROUTINE HPCG DOES NOT REQUIRE
C                    AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL
C                    FOR HANDING RESULT VALUES TO THE MAIN PROGRAM
C                    (CALLING HPCG) WHICH ARE OBTAINED BY SPECIAL
C                    MANIPULATIONS WITH OUTPUT DATA IN SUBROUTINE OUTP.
C           Y      - INPUT VECTOR OF INITIAL VALUES.  (DESTROYED)
C                    LATERON Y IS THE RESULTING VECTOR OF DEPENDENT
C                    VARIABLES COMPUTED AT INTERMEDIATE POINTS X.
C           DERY   - INPUT VECTOR OF ERROR WEIGHTS.  (DESTROYED)
C                    THE SUM OF ITS COMPONENTS MUST BE EQUAL TO 1.
C                    LATERON DERY IS THE VECTOR OF DERIVATIVES, WHICH
C                    BELONG TO FUNCTION VALUES Y AT A POINT X.
C           NDIM   - AN INPUT VALUE, WHICH SPECIFIES THE NUMBER OF
C                    EQUATIONS IN THE SYSTEM.
C           IHLF   - AN OUTPUT VALUE, WHICH SPECIFIES THE NUMBER OF
C                    BISECTIONS OF THE INITIAL INCREMENT. IF IHLF GETS
C                    GREATER THAN 10, SUBROUTINE HPCG RETURNS WITH
C                    ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM.
C                    ERROR MESSAGE IHLF=12 OR IHLF=13 APPEARS IN CASE
C                    PRMT(3)=0 OR IN CASE SIGN(PRMT(3)).NE.SIGN(PRMT(2)-
C                    PRMT(1)) RESPECTIVELY.
C           FCT    - THE NAME OF AN EXTERNAL SUBROUTINE USED. IT
C                    COMPUTES THE RIGHT HAND SIDES DERY OF THE SYSTEM
C                    TO GIVEN VALUES OF X AND Y. ITS PARAMETER LIST
C                    MUST BE X,Y,DERY. THE SUBROUTINE SHOULD NOT
C                    DESTROY X AND Y.
C           OUTP   - THE NAME OF AN EXTERNAL OUTPUT SUBROUTINE USED.
C                    ITS PARAMETER LIST MUST BE X,Y,DERY,IHLF,NDIM,PRMT.
C                    NONE OF THESE PARAMETERS (EXCEPT, IF NECESSARY,
C                    PRMT(4),PRMT(5),...) SHOULD BE CHANGED BY
C                    SUBROUTINE OUTP. IF PRMT(5) IS CHANGED TO NON-ZERO,
C                    SUBROUTINE HPCG IS TERMINATED.
C           AUX    - AN AUXILIARY STORAGE ARRAY WITH 16 ROWS AND NDIM
C                    COLUMNS.
C
C        REMARKS
C           THE PROCEDURE TERMINATES AND RETURNS TO CALLING PROGRAM, IF
C           (1) MORE THAN 10 BISECTIONS OF THE INITIAL INCREMENT ARE
C               NECESSARY TO GET SATISFACTORY ACCURACY (ERROR MESSAGE
C               IHLF=11),
C           (2) INITIAL INCREMENT IS EQUAL TO 0 OR HAS WRONG SIGN
C               (ERROR MESSAGES IHLF=12 OR IHLF=13),
C           (3) THE WHOLE INTEGRATION INTERVAL IS WORKED THROUGH,
C           (4) SUBROUTINE OUTP HAS CHANGED PRMT(5) TO NON-ZERO.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           THE EXTERNAL SUBROUTINES FCT(X,Y,DERY) AND
C           OUTP(X,Y,DERY,IHLF,NDIM,PRMT) MUST BE FURNISHED BY THE USER.
C
C        METHOD
C           EVALUATION IS DONE BY MEANS OF HAMMINGS MODIFIED PREDICTOR-
C           CORRECTOR METHOD. IT IS A FOURTH ORDER METHOD, USING 4
C           PRECEEDING POINTS FOR COMPUTATION OF A NEW VECTOR Y OF THE
C           DEPENDENT VARIABLES.
C           FOURTH ORDER RUNGE-KUTTA METHOD SUGGESTED BY RALSTON IS
C           USED FOR ADJUSTMENT OF THE INITIAL INCREMENT AND FOR
C           COMPUTATION OF STARTING VALUES.
C           SUBROUTINE HPCG AUTOMATICALLY ADJUSTS THE INCREMENT DURING
C           THE WHOLE COMPUTATION BY HALVING OR DOUBLING.
C           TO GET FULL FLEXIBILITY IN OUTPUT, AN OUTPUT SUBROUTINE
C           MUST BE CODED BY THE USER.
C           FOR REFERENCE, SEE
C           (1)  RALSTON/WILF, MATHEMATICAL METHODS FOR DIGITAL
C                COMPUTERS, WILEY, NEW YORK/LONDON, 1960, PP.95-109.
C           (2)  RALSTON, RUNGE-KUTTA METHODS WITH MINIMUM ERROR BOUNDS,
C                MTAC, VOL.16, ISS.80 (1962), PP.431-437.
C
C     ..................................................................
C
      SUBROUTINE HPCG(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)
C
      DIMENSION PRMT(*),Y(*),DERY(*),AUX(16,*)
      N=1
      IHLF=0
      X=PRMT(1)
      H=PRMT(3)
      PRMT(5)=0.
      DO 1 I=1,NDIM
      AUX(16,I)=0.
      AUX(15,I)=DERY(I)
    1 AUX(1,I)=Y(I)
      IF(H*(PRMT(2)-X))3,2,4
C
C     ERROR RETURNS
    2 IHLF=12
      GOTO 4
    3 IHLF=13
C
C     COMPUTATION OF DERY FOR STARTING VALUES
    4 CALL FCT(X,Y,DERY)
C
C     RECORDING OF STARTING VALUES
      CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)
      IF(PRMT(5))6,5,6
    5 IF(IHLF)7,7,6
    6 RETURN
    7 DO 8 I=1,NDIM
    8 AUX(8,I)=DERY(I)
C
C     COMPUTATION OF AUX(2,I)
      ISW=1
      GOTO 100
C
    9 X=X+H
      DO 10 I=1,NDIM
   10 AUX(2,I)=Y(I)
C
C     INCREMENT H IS TESTED BY MEANS OF BISECTION
   11 IHLF=IHLF+1
      X=X-H
      DO 12 I=1,NDIM
   12 AUX(4,I)=AUX(2,I)
      H=.5*H
      N=1
      ISW=2
      GOTO 100
C
   13 X=X+H
      CALL FCT(X,Y,DERY)
      N=2
      DO 14 I=1,NDIM
      AUX(2,I)=Y(I)
   14 AUX(9,I)=DERY(I)
      ISW=3
      GOTO 100
C
C     COMPUTATION OF TEST VALUE DELT
   15 DELT=0.
      DO 16 I=1,NDIM
   16 DELT=DELT+AUX(15,I)*ABS(Y(I)-AUX(4,I))
      DELT=.06666667*DELT
      IF(DELT-PRMT(4))19,19,17
   17 IF(IHLF-10)11,18,18
C
C     NO SATISFACTORY ACCURACY AFTER 10 BISECTIONS. ERROR MESSAGE.
   18 IHLF=11
      X=X+H
      GOTO 4
C
C     THERE IS SATISFACTORY ACCURACY AFTER LESS THAN 11 BISECTIONS.
   19 X=X+H
      CALL FCT(X,Y,DERY)
      DO 20 I=1,NDIM
      AUX(3,I)=Y(I)
   20 AUX(10,I)=DERY(I)
      N=3
      ISW=4
      GOTO 100
C
   21 N=1
      X=X+H
      CALL FCT(X,Y,DERY)
      X=PRMT(1)
      DO 22 I=1,NDIM
      AUX(11,I)=DERY(I)
   220Y(I)=AUX(1,I)+H*(.375*AUX(8,I)+.7916667*AUX(9,I)
     1-.2083333*AUX(10,I)+.04166667*DERY(I))
   23 X=X+H
      N=N+1
      CALL FCT(X,Y,DERY)
      CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)
      IF(PRMT(5))6,24,6
   24 IF(N-4)25,200,200
   25 DO 26 I=1,NDIM
      AUX(N,I)=Y(I)
   26 AUX(N+7,I)=DERY(I)
      IF(N-3)27,29,200
C
   27 DO 28 I=1,NDIM
      DELT=AUX(9,I)+AUX(9,I)
      DELT=DELT+DELT
   28 Y(I)=AUX(1,I)+.3333333*H*(AUX(8,I)+DELT+AUX(10,I))
      GOTO 23
C
   29 DO 30 I=1,NDIM
      DELT=AUX(9,I)+AUX(10,I)
      DELT=DELT+DELT+DELT
   30 Y(I)=AUX(1,I)+.375*H*(AUX(8,I)+DELT+AUX(11,I))
      GOTO 23
C
C     THE FOLLOWING PART OF SUBROUTINE HPCG COMPUTES BY MEANS OF
C     RUNGE-KUTTA METHOD STARTING VALUES FOR THE NOT SELF-STARTING
C     PREDICTOR-CORRECTOR METHOD.
  100 DO 101 I=1,NDIM
      Z=H*AUX(N+7,I)
      AUX(5,I)=Z
  101 Y(I)=AUX(N,I)+.4*Z
C     Z IS AN AUXILIARY STORAGE LOCATION
C
      Z=X+.4*H
      CALL FCT(Z,Y,DERY)
      DO 102 I=1,NDIM
      Z=H*DERY(I)
      AUX(6,I)=Z
  102 Y(I)=AUX(N,I)+.2969776*AUX(5,I)+.1587596*Z
C
      Z=X+.4557372*H
      CALL FCT(Z,Y,DERY)
      DO 103 I=1,NDIM
      Z=H*DERY(I)
      AUX(7,I)=Z
  103 Y(I)=AUX(N,I)+.2181004*AUX(5,I)-3.050965*AUX(6,I)+3.832865*Z
C
      Z=X+H
      CALL FCT(Z,Y,DERY)
      DO 104 I=1,NDIM
  1040Y(I)=AUX(N,I)+.1747603*AUX(5,I)-.5514807*AUX(6,I)
     1+1.205536*AUX(7,I)+.1711848*H*DERY(I)
      GOTO(9,13,15,21),ISW
C
C     POSSIBLE BREAK-POINT FOR LINKAGE
C
C     STARTING VALUES ARE COMPUTED.
C     NOW START HAMMINGS MODIFIED PREDICTOR-CORRECTOR METHOD.
  200 ISTEP=3
  201 IF(N-8)204,202,204
C
C     N=8 CAUSES THE ROWS OF AUX TO CHANGE THEIR STORAGE LOCATIONS
  202 DO 203 N=2,7
      DO 203 I=1,NDIM
      AUX(N-1,I)=AUX(N,I)
  203 AUX(N+6,I)=AUX(N+7,I)
      N=7
C
C     N LESS THAN 8 CAUSES N+1 TO GET N
  204 N=N+1
C
C     COMPUTATION OF NEXT VECTOR Y
      DO 205 I=1,NDIM
      AUX(N-1,I)=Y(I)
  205 AUX(N+6,I)=DERY(I)
      X=X+H
  206 ISTEP=ISTEP+1
      DO 207 I=1,NDIM
     0DELT=AUX(N-4,I)+1.333333*H*(AUX(N+6,I)+AUX(N+6,I)-AUX(N+5,I)+
     1AUX(N+4,I)+AUX(N+4,I))
      Y(I)=DELT-.9256198*AUX(16,I)
  207 AUX(16,I)=DELT
C     PREDICTOR IS NOW GENERATED IN ROW 16 OF AUX, MODIFIED PREDICTOR
C     IS GENERATED IN Y. DELT MEANS AN AUXILIARY STORAGE.
C
      CALL FCT(X,Y,DERY)
C     DERIVATIVE OF MODIFIED PREDICTOR IS GENERATED IN DERY
C
      DO 208 I=1,NDIM
     0DELT=.125*(9.*AUX(N-1,I)-AUX(N-3,I)+3.*H*(DERY(I)+AUX(N+6,I)+
     1AUX(N+6,I)-AUX(N+5,I)))
      AUX(16,I)=AUX(16,I)-DELT
  208 Y(I)=DELT+.07438017*AUX(16,I)
C
C     TEST WHETHER H MUST BE HALVED OR DOUBLED
      DELT=0.
      DO 209 I=1,NDIM
  209 DELT=DELT+AUX(15,I)*ABS(AUX(16,I))
      IF(DELT-PRMT(4))210,222,222
C
C     H MUST NOT BE HALVED. THAT MEANS Y(I) ARE GOOD.
  210 CALL FCT(X,Y,DERY)
      CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)
      IF(PRMT(5))212,211,212
  211 IF(IHLF-11)213,212,212
  212 RETURN
  213 IF(H*(X-PRMT(2)))214,212,212
  214 IF(ABS(X-PRMT(2))-.1*ABS(H))212,215,215
  215 IF(DELT-.02*PRMT(4))216,216,201
C
C
C     H COULD BE DOUBLED IF ALL NECESSARY PRECEEDING VALUES ARE
C     AVAILABLE
  216 IF(IHLF)201,201,217
  217 IF(N-7)201,218,218
  218 IF(ISTEP-4)201,219,219
  219 IMOD=ISTEP/2
      IF(ISTEP-IMOD-IMOD)201,220,201
  220 H=H+H
      IHLF=IHLF-1
      ISTEP=0
      DO 221 I=1,NDIM
      AUX(N-1,I)=AUX(N-2,I)
      AUX(N-2,I)=AUX(N-4,I)
      AUX(N-3,I)=AUX(N-6,I)
      AUX(N+6,I)=AUX(N+5,I)
      AUX(N+5,I)=AUX(N+3,I)
      AUX(N+4,I)=AUX(N+1,I)
      DELT=AUX(N+6,I)+AUX(N+5,I)
      DELT=DELT+DELT+DELT
  2210AUX(16,I)=8.962963*(Y(I)-AUX(N-3,I))-3.361111*H*(DERY(I)+DELT
     1+AUX(N+4,I))
      GOTO 201
C
C
C     H MUST BE HALVED
  222 IHLF=IHLF+1
      IF(IHLF-10)223,223,210
  223 H=.5*H
      ISTEP=0
      DO 224 I=1,NDIM
     0Y(I)=.00390625*(80.*AUX(N-1,I)+135.*AUX(N-2,I)+40.*AUX(N-3,I)+
     1AUX(N-4,I))-.1171875*(AUX(N+6,I)-6.*AUX(N+5,I)-AUX(N+4,I))*H
     0AUX(N-4,I)=.00390625*(12.*AUX(N-1,I)+135.*AUX(N-2,I)+
     1108.*AUX(N-3,I)+AUX(N-4,I))-.0234375*(AUX(N+6,I)+18.*AUX(N+5,I)-
     29.*AUX(N+4,I))*H
      AUX(N-3,I)=AUX(N-2,I)
  224 AUX(N+4,I)=AUX(N+5,I)
      X=X-H
      DELT=X-(H+H)
      CALL FCT(DELT,Y,DERY)
      DO 225 I=1,NDIM
      AUX(N-2,I)=Y(I)
      AUX(N+5,I)=DERY(I)
  225 Y(I)=AUX(N-4,I)
      DELT=DELT-(H+H)
      CALL FCT(DELT,Y,DERY)
      DO 226 I=1,NDIM
      DELT=AUX(N+5,I)+AUX(N+4,I)
      DELT=DELT+DELT+DELT
     0AUX(16,I)=8.962963*(AUX(N-1,I)-Y(I))-3.361111*H*(AUX(N+6,I)+DELT
     1+DERY(I))
  226 AUX(N+3,I)=DERY(I)
      GOTO 206
      END
C
C=======================================================================
C