Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/hpcg.ssp
There are 2 other files named hpcg.ssp in the archive. Click here to see a list.
C                                                                       HPCG  10
C                                                                       HPCG  20
C     ..................................................................HPCG  30
C                                                                       HPCG  40
C        SUBROUTINE HPCG                                                HPCG  50
C                                                                       HPCG  60
C        PURPOSE                                                        HPCG  70
C           TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY GENERAL           HPCG  80
C           DIFFERENTIAL EQUATIONS WITH GIVEN INITIAL VALUES.           HPCG  90
C                                                                       HPCG 100
C        USAGE                                                          HPCG 110
C           CALL HPCG (PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)              HPCG 120
C           PARAMETERS FCT AND OUTP REQUIRE AN EXTERNAL STATEMENT.      HPCG 130
C                                                                       HPCG 140
C        DESCRIPTION OF PARAMETERS                                      HPCG 150
C           PRMT   - AN INPUT AND OUTPUT VECTOR WITH DIMENSION GREATER  HPCG 160
C                    OR EQUAL TO 5, WHICH SPECIFIES THE PARAMETERS OF   HPCG 170
C                    THE INTERVAL AND OF ACCURACY AND WHICH SERVES FOR  HPCG 180
C                    COMMUNICATION BETWEEN OUTPUT SUBROUTINE (FURNISHED HPCG 190
C                    BY THE USER) AND SUBROUTINE HPCG. EXCEPT PRMT(5)   HPCG 200
C                    THE COMPONENTS ARE NOT DESTROYED BY SUBROUTINE     HPCG 210
C                    HPCG AND THEY ARE                                  HPCG 220
C           PRMT(1)- LOWER BOUND OF THE INTERVAL (INPUT),               HPCG 230
C           PRMT(2)- UPPER BOUND OF THE INTERVAL (INPUT),               HPCG 240
C           PRMT(3)- INITIAL INCREMENT OF THE INDEPENDENT VARIABLE      HPCG 250
C                    (INPUT),                                           HPCG 260
C           PRMT(4)- UPPER ERROR BOUND (INPUT). IF ABSOLUTE ERROR IS    HPCG 270
C                    GREATER THAN PRMT(4), INCREMENT GETS HALVED.       HPCG 280
C                    IF INCREMENT IS LESS THAN PRMT(3) AND ABSOLUTE     HPCG 290
C                    ERROR LESS THAN PRMT(4)/50, INCREMENT GETS DOUBLED.HPCG 300
C                    THE USER MAY CHANGE PRMT(4) BY MEANS OF HIS        HPCG 310
C                    OUTPUT SUBROUTINE.                                 HPCG 320
C           PRMT(5)- NO INPUT PARAMETER. SUBROUTINE HPCG INITIALIZES    HPCG 330
C                    PRMT(5)=0. IF THE USER WANTS TO TERMINATE          HPCG 340
C                    SUBROUTINE HPCG AT ANY OUTPUT POINT, HE HAS TO     HPCG 350
C                    CHANGE PRMT(5) TO NON-ZERO BY MEANS OF SUBROUTINE  HPCG 360
C                    OUTP. FURTHER COMPONENTS OF VECTOR PRMT ARE        HPCG 370
C                    FEASIBLE IF ITS DIMENSION IS DEFINED GREATER       HPCG 380
C                    THAN 5. HOWEVER SUBROUTINE HPCG DOES NOT REQUIRE   HPCG 390
C                    AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL   HPCG 400
C                    FOR HANDING RESULT VALUES TO THE MAIN PROGRAM      HPCG 410
C                    (CALLING HPCG) WHICH ARE OBTAINED BY SPECIAL       HPCG 420
C                    MANIPULATIONS WITH OUTPUT DATA IN SUBROUTINE OUTP. HPCG 430
C           Y      - INPUT VECTOR OF INITIAL VALUES.  (DESTROYED)       HPCG 440
C                    LATERON Y IS THE RESULTING VECTOR OF DEPENDENT     HPCG 450
C                    VARIABLES COMPUTED AT INTERMEDIATE POINTS X.       HPCG 460
C           DERY   - INPUT VECTOR OF ERROR WEIGHTS.  (DESTROYED)        HPCG 470
C                    THE SUM OF ITS COMPONENTS MUST BE EQUAL TO 1.      HPCG 480
C                    LATERON DERY IS THE VECTOR OF DERIVATIVES, WHICH   HPCG 490
C                    BELONG TO FUNCTION VALUES Y AT A POINT X.          HPCG 500
C           NDIM   - AN INPUT VALUE, WHICH SPECIFIES THE NUMBER OF      HPCG 510
C                    EQUATIONS IN THE SYSTEM.                           HPCG 520
C           IHLF   - AN OUTPUT VALUE, WHICH SPECIFIES THE NUMBER OF     HPCG 530
C                    BISECTIONS OF THE INITIAL INCREMENT. IF IHLF GETS  HPCG 540
C                    GREATER THAN 10, SUBROUTINE HPCG RETURNS WITH      HPCG 550
C                    ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM.           HPCG 560
C                    ERROR MESSAGE IHLF=12 OR IHLF=13 APPEARS IN CASE   HPCG 570
C                    PRMT(3)=0 OR IN CASE SIGN(PRMT(3)).NE.SIGN(PRMT(2)-HPCG 580
C                    PRMT(1)) RESPECTIVELY.                             HPCG 590
C           FCT    - THE NAME OF AN EXTERNAL SUBROUTINE USED. IT        HPCG 600
C                    COMPUTES THE RIGHT HAND SIDES DERY OF THE SYSTEM   HPCG 610
C                    TO GIVEN VALUES OF X AND Y. ITS PARAMETER LIST     HPCG 620
C                    MUST BE X,Y,DERY. THE SUBROUTINE SHOULD NOT        HPCG 630
C                    DESTROY X AND Y.                                   HPCG 640
C           OUTP   - THE NAME OF AN EXTERNAL OUTPUT SUBROUTINE USED.    HPCG 650
C                    ITS PARAMETER LIST MUST BE X,Y,DERY,IHLF,NDIM,PRMT.HPCG 660
C                    NONE OF THESE PARAMETERS (EXCEPT, IF NECESSARY,    HPCG 670
C                    PRMT(4),PRMT(5),...) SHOULD BE CHANGED BY          HPCG 680
C                    SUBROUTINE OUTP. IF PRMT(5) IS CHANGED TO NON-ZERO,HPCG 690
C                    SUBROUTINE HPCG IS TERMINATED.                     HPCG 700
C           AUX    - AN AUXILIARY STORAGE ARRAY WITH 16 ROWS AND NDIM   HPCG 710
C                    COLUMNS.                                           HPCG 720
C                                                                       HPCG 730
C        REMARKS                                                        HPCG 740
C           THE PROCEDURE TERMINATES AND RETURNS TO CALLING PROGRAM, IF HPCG 750
C           (1) MORE THAN 10 BISECTIONS OF THE INITIAL INCREMENT ARE    HPCG 760
C               NECESSARY TO GET SATISFACTORY ACCURACY (ERROR MESSAGE   HPCG 770
C               IHLF=11),                                               HPCG 780
C           (2) INITIAL INCREMENT IS EQUAL TO 0 OR HAS WRONG SIGN       HPCG 790
C               (ERROR MESSAGES IHLF=12 OR IHLF=13),                    HPCG 800
C           (3) THE WHOLE INTEGRATION INTERVAL IS WORKED THROUGH,       HPCG 810
C           (4) SUBROUTINE OUTP HAS CHANGED PRMT(5) TO NON-ZERO.        HPCG 820
C                                                                       HPCG 830
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  HPCG 840
C           THE EXTERNAL SUBROUTINES FCT(X,Y,DERY) AND                  HPCG 850
C           OUTP(X,Y,DERY,IHLF,NDIM,PRMT) MUST BE FURNISHED BY THE USER.HPCG 860
C                                                                       HPCG 870
C        METHOD                                                         HPCG 880
C           EVALUATION IS DONE BY MEANS OF HAMMINGS MODIFIED PREDICTOR- HPCG 890
C           CORRECTOR METHOD. IT IS A FOURTH ORDER METHOD, USING 4      HPCG 900
C           PRECEEDING POINTS FOR COMPUTATION OF A NEW VECTOR Y OF THE  HPCG 910
C           DEPENDENT VARIABLES.                                        HPCG 920
C           FOURTH ORDER RUNGE-KUTTA METHOD SUGGESTED BY RALSTON IS     HPCG 930
C           USED FOR ADJUSTMENT OF THE INITIAL INCREMENT AND FOR        HPCG 940
C           COMPUTATION OF STARTING VALUES.                             HPCG 950
C           SUBROUTINE HPCG AUTOMATICALLY ADJUSTS THE INCREMENT DURING  HPCG 960
C           THE WHOLE COMPUTATION BY HALVING OR DOUBLING.               HPCG 970
C           TO GET FULL FLEXIBILITY IN OUTPUT, AN OUTPUT SUBROUTINE     HPCG 980
C           MUST BE CODED BY THE USER.                                  HPCG 990
C           FOR REFERENCE, SEE                                          HPCG1000
C           (1)  RALSTON/WILF, MATHEMATICAL METHODS FOR DIGITAL         HPCG1010
C                COMPUTERS, WILEY, NEW YORK/LONDON, 1960, PP.95-109.    HPCG1020
C           (2)  RALSTON, RUNGE-KUTTA METHODS WITH MINIMUM ERROR BOUNDS,HPCG1030
C                MTAC, VOL.16, ISS.80 (1962), PP.431-437.               HPCG1040
C                                                                       HPCG1050
C     ..................................................................HPCG1060
C                                                                       HPCG1070
      SUBROUTINE HPCG(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               HPCG1080
C                                                                       HPCG1090
C                                                                       HPCG1100
      DIMENSION PRMT(1),Y(1),DERY(1),AUX(16,1)                          HPCG1110
      N=1                                                               HPCG1120
      IHLF=0                                                            HPCG1130
      X=PRMT(1)                                                         HPCG1140
      H=PRMT(3)                                                         HPCG1150
      PRMT(5)=0.                                                        HPCG1160
      DO 1 I=1,NDIM                                                     HPCG1170
      AUX(16,I)=0.                                                      HPCG1180
      AUX(15,I)=DERY(I)                                                 HPCG1190
    1 AUX(1,I)=Y(I)                                                     HPCG1200
      IF(H*(PRMT(2)-X))3,2,4                                            HPCG1210
C                                                                       HPCG1220
C     ERROR RETURNS                                                     HPCG1230
    2 IHLF=12                                                           HPCG1240
      GOTO 4                                                            HPCG1250
    3 IHLF=13                                                           HPCG1260
C                                                                       HPCG1270
C     COMPUTATION OF DERY FOR STARTING VALUES                           HPCG1280
    4 CALL FCT(X,Y,DERY)                                                HPCG1290
C                                                                       HPCG1300
C     RECORDING OF STARTING VALUES                                      HPCG1310
      CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                HPCG1320
      IF(PRMT(5))6,5,6                                                  HPCG1330
    5 IF(IHLF)7,7,6                                                     HPCG1340
    6 RETURN                                                            HPCG1350
    7 DO 8 I=1,NDIM                                                     HPCG1360
    8 AUX(8,I)=DERY(I)                                                  HPCG1370
C                                                                       HPCG1380
C     COMPUTATION OF AUX(2,I)                                           HPCG1390
      ISW=1                                                             HPCG1400
      GOTO 100                                                          HPCG1410
C                                                                       HPCG1420
    9 X=X+H                                                             HPCG1430
      DO 10 I=1,NDIM                                                    HPCG1440
   10 AUX(2,I)=Y(I)                                                     HPCG1450
C                                                                       HPCG1460
C     INCREMENT H IS TESTED BY MEANS OF BISECTION                       HPCG1470
   11 IHLF=IHLF+1                                                       HPCG1480
      X=X-H                                                             HPCG1490
      DO 12 I=1,NDIM                                                    HPCG1500
   12 AUX(4,I)=AUX(2,I)                                                 HPCG1510
      H=.5*H                                                            HPCG1520
      N=1                                                               HPCG1530
      ISW=2                                                             HPCG1540
      GOTO 100                                                          HPCG1550
C                                                                       HPCG1560
   13 X=X+H                                                             HPCG1570
      CALL FCT(X,Y,DERY)                                                HPCG1580
      N=2                                                               HPCG1590
      DO 14 I=1,NDIM                                                    HPCG1600
      AUX(2,I)=Y(I)                                                     HPCG1610
   14 AUX(9,I)=DERY(I)                                                  HPCG1620
      ISW=3                                                             HPCG1630
      GOTO 100                                                          HPCG1640
C                                                                       HPCG1650
C     COMPUTATION OF TEST VALUE DELT                                    HPCG1660
   15 DELT=0.                                                           HPCG1670
      DO 16 I=1,NDIM                                                    HPCG1680
   16 DELT=DELT+AUX(15,I)*ABS(Y(I)-AUX(4,I))                            HPCG1690
      DELT=.06666667*DELT                                               HPCG1700
      IF(DELT-PRMT(4))19,19,17                                          HPCG1710
   17 IF(IHLF-10)11,18,18                                               HPCG1720
C                                                                       HPCG1730
C     NO SATISFACTORY ACCURACY AFTER 10 BISECTIONS. ERROR MESSAGE.      HPCG1740
   18 IHLF=11                                                           HPCG1750
      X=X+H                                                             HPCG1760
      GOTO 4                                                            HPCG1770
C                                                                       HPCG1780
C     THERE IS SATISFACTORY ACCURACY AFTER LESS THAN 11 BISECTIONS.     HPCG1790
   19 X=X+H                                                             HPCG1800
      CALL FCT(X,Y,DERY)                                                HPCG1810
      DO 20 I=1,NDIM                                                    HPCG1820
      AUX(3,I)=Y(I)                                                     HPCG1830
   20 AUX(10,I)=DERY(I)                                                 HPCG1840
      N=3                                                               HPCG1850
      ISW=4                                                             HPCG1860
      GOTO 100                                                          HPCG1870
C                                                                       HPCG1880
   21 N=1                                                               HPCG1890
      X=X+H                                                             HPCG1900
      CALL FCT(X,Y,DERY)                                                HPCG1910
      X=PRMT(1)                                                         HPCG1920
      DO 22 I=1,NDIM                                                    HPCG1930
      AUX(11,I)=DERY(I)                                                 HPCG1940
   220Y(I)=AUX(1,I)+H*(.375*AUX(8,I)+.7916667*AUX(9,I)                  HPCG1950
     1-.2083333*AUX(10,I)+.04166667*DERY(I))                            HPCG1960
   23 X=X+H                                                             HPCG1970
      N=N+1                                                             HPCG1980
      CALL FCT(X,Y,DERY)                                                HPCG1990
      CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                HPCG2000
      IF(PRMT(5))6,24,6                                                 HPCG2010
   24 IF(N-4)25,200,200                                                 HPCG2020
   25 DO 26 I=1,NDIM                                                    HPCG2030
      AUX(N,I)=Y(I)                                                     HPCG2040
   26 AUX(N+7,I)=DERY(I)                                                HPCG2050
      IF(N-3)27,29,200                                                  HPCG2060
C                                                                       HPCG2070
   27 DO 28 I=1,NDIM                                                    HPCG2080
      DELT=AUX(9,I)+AUX(9,I)                                            HPCG2090
      DELT=DELT+DELT                                                    HPCG2100
   28 Y(I)=AUX(1,I)+.3333333*H*(AUX(8,I)+DELT+AUX(10,I))                HPCG2110
      GOTO 23                                                           HPCG2120
C                                                                       HPCG2130
   29 DO 30 I=1,NDIM                                                    HPCG2140
      DELT=AUX(9,I)+AUX(10,I)                                           HPCG2150
      DELT=DELT+DELT+DELT                                               HPCG2160
   30 Y(I)=AUX(1,I)+.375*H*(AUX(8,I)+DELT+AUX(11,I))                    HPCG2170
      GOTO 23                                                           HPCG2180
C                                                                       HPCG2190
C     THE FOLLOWING PART OF SUBROUTINE HPCG COMPUTES BY MEANS OF        HPCG2200
C     RUNGE-KUTTA METHOD STARTING VALUES FOR THE NOT SELF-STARTING      HPCG2210
C     PREDICTOR-CORRECTOR METHOD.                                       HPCG2220
  100 DO 101 I=1,NDIM                                                   HPCG2230
      Z=H*AUX(N+7,I)                                                    HPCG2240
      AUX(5,I)=Z                                                        HPCG2250
  101 Y(I)=AUX(N,I)+.4*Z                                                HPCG2260
C     Z IS AN AUXILIARY STORAGE LOCATION                                HPCG2270
C                                                                       HPCG2280
      Z=X+.4*H                                                          HPCG2290
      CALL FCT(Z,Y,DERY)                                                HPCG2300
      DO 102 I=1,NDIM                                                   HPCG2310
      Z=H*DERY(I)                                                       HPCG2320
      AUX(6,I)=Z                                                        HPCG2330
  102 Y(I)=AUX(N,I)+.2969776*AUX(5,I)+.1587596*Z                        HPCG2340
C                                                                       HPCG2350
      Z=X+.4557372*H                                                    HPCG2360
      CALL FCT(Z,Y,DERY)                                                HPCG2370
      DO 103 I=1,NDIM                                                   HPCG2380
      Z=H*DERY(I)                                                       HPCG2390
      AUX(7,I)=Z                                                        HPCG2400
  103 Y(I)=AUX(N,I)+.2181004*AUX(5,I)-3.050965*AUX(6,I)+3.832865*Z      HPCG2410
C                                                                       HPCG2420
      Z=X+H                                                             HPCG2430
      CALL FCT(Z,Y,DERY)                                                HPCG2440
      DO 104 I=1,NDIM                                                   HPCG2450
  1040Y(I)=AUX(N,I)+.1747603*AUX(5,I)-.5514807*AUX(6,I)                 HPCG2460
     1+1.205536*AUX(7,I)+.1711848*H*DERY(I)                             HPCG2470
      GOTO(9,13,15,21),ISW                                              HPCG2480
C                                                                       HPCG2490
C     POSSIBLE BREAK-POINT FOR LINKAGE                                  HPCG2500
C                                                                       HPCG2510
C     STARTING VALUES ARE COMPUTED.                                     HPCG2520
C     NOW START HAMMINGS MODIFIED PREDICTOR-CORRECTOR METHOD.           HPCG2530
  200 ISTEP=3                                                           HPCG2540
  201 IF(N-8)204,202,204                                                HPCG2550
C                                                                       HPCG2560
C     N=8 CAUSES THE ROWS OF AUX TO CHANGE THEIR STORAGE LOCATIONS      HPCG2570
  202 DO 203 N=2,7                                                      HPCG2580
      DO 203 I=1,NDIM                                                   HPCG2590
      AUX(N-1,I)=AUX(N,I)                                               HPCG2600
  203 AUX(N+6,I)=AUX(N+7,I)                                             HPCG2610
      N=7                                                               HPCG2620
C                                                                       HPCG2630
C     N LESS THAN 8 CAUSES N+1 TO GET N                                 HPCG2640
  204 N=N+1                                                             HPCG2650
C                                                                       HPCG2660
C     COMPUTATION OF NEXT VECTOR Y                                      HPCG2670
      DO 205 I=1,NDIM                                                   HPCG2680
      AUX(N-1,I)=Y(I)                                                   HPCG2690
  205 AUX(N+6,I)=DERY(I)                                                HPCG2700
      X=X+H                                                             HPCG2710
  206 ISTEP=ISTEP+1                                                     HPCG2720
      DO 207 I=1,NDIM                                                   HPCG2730
     0DELT=AUX(N-4,I)+1.333333*H*(AUX(N+6,I)+AUX(N+6,I)-AUX(N+5,I)+     HPCG2740
     1AUX(N+4,I)+AUX(N+4,I))                                            HPCG2750
      Y(I)=DELT-.9256198*AUX(16,I)                                      HPCG2760
  207 AUX(16,I)=DELT                                                    HPCG2770
C     PREDICTOR IS NOW GENERATED IN ROW 16 OF AUX, MODIFIED PREDICTOR   HPCG2780
C     IS GENERATED IN Y. DELT MEANS AN AUXILIARY STORAGE.               HPCG2790
C                                                                       HPCG2800
      CALL FCT(X,Y,DERY)                                                HPCG2810
C     DERIVATIVE OF MODIFIED PREDICTOR IS GENERATED IN DERY             HPCG2820
C                                                                       HPCG2830
      DO 208 I=1,NDIM                                                   HPCG2840
     0DELT=.125*(9.*AUX(N-1,I)-AUX(N-3,I)+3.*H*(DERY(I)+AUX(N+6,I)+     HPCG2850
     1AUX(N+6,I)-AUX(N+5,I)))                                           HPCG2860
      AUX(16,I)=AUX(16,I)-DELT                                          HPCG2870
  208 Y(I)=DELT+.07438017*AUX(16,I)                                     HPCG2880
C                                                                       HPCG2890
C     TEST WHETHER H MUST BE HALVED OR DOUBLED                          HPCG2900
      DELT=0.                                                           HPCG2910
      DO 209 I=1,NDIM                                                   HPCG2920
  209 DELT=DELT+AUX(15,I)*ABS(AUX(16,I))                                HPCG2930
      IF(DELT-PRMT(4))210,222,222                                       HPCG2940
C                                                                       HPCG2950
C     H MUST NOT BE HALVED. THAT MEANS Y(I) ARE GOOD.                   HPCG2960
  210 CALL FCT(X,Y,DERY)                                                HPCG2970
      CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                HPCG2980
      IF(PRMT(5))212,211,212                                            HPCG2990
  211 IF(IHLF-11)213,212,212                                            HPCG3000
  212 RETURN                                                            HPCG3010
  213 IF(H*(X-PRMT(2)))214,212,212                                      HPCG3020
  214 IF(ABS(X-PRMT(2))-.1*ABS(H))212,215,215                           HPCG3030
  215 IF(DELT-.02*PRMT(4))216,216,201                                   HPCG3040
C                                                                       HPCG3050
C                                                                       HPCG3060
C     H COULD BE DOUBLED IF ALL NECESSARY PRECEEDING VALUES ARE         HPCG3070
C     AVAILABLE                                                         HPCG3080
  216 IF(IHLF)201,201,217                                               HPCG3090
  217 IF(N-7)201,218,218                                                HPCG3100
  218 IF(ISTEP-4)201,219,219                                            HPCG3110
  219 IMOD=ISTEP/2                                                      HPCG3120
      IF(ISTEP-IMOD-IMOD)201,220,201                                    HPCG3130
  220 H=H+H                                                             HPCG3140
      IHLF=IHLF-1                                                       HPCG3150
      ISTEP=0                                                           HPCG3160
      DO 221 I=1,NDIM                                                   HPCG3170
      AUX(N-1,I)=AUX(N-2,I)                                             HPCG3180
      AUX(N-2,I)=AUX(N-4,I)                                             HPCG3190
      AUX(N-3,I)=AUX(N-6,I)                                             HPCG3200
      AUX(N+6,I)=AUX(N+5,I)                                             HPCG3210
      AUX(N+5,I)=AUX(N+3,I)                                             HPCG3220
      AUX(N+4,I)=AUX(N+1,I)                                             HPCG3230
      DELT=AUX(N+6,I)+AUX(N+5,I)                                        HPCG3240
      DELT=DELT+DELT+DELT                                               HPCG3250
  2210AUX(16,I)=8.962963*(Y(I)-AUX(N-3,I))-3.361111*H*(DERY(I)+DELT     HPCG3260
     1+AUX(N+4,I))                                                      HPCG3270
      GOTO 201                                                          HPCG3280
C                                                                       HPCG3290
C                                                                       HPCG3300
C     H MUST BE HALVED                                                  HPCG3310
  222 IHLF=IHLF+1                                                       HPCG3320
      IF(IHLF-10)223,223,210                                            HPCG3330
  223 H=.5*H                                                            HPCG3340
      ISTEP=0                                                           HPCG3350
      DO 224 I=1,NDIM                                                   HPCG3360
     0Y(I)=.00390625*(80.*AUX(N-1,I)+135.*AUX(N-2,I)+40.*AUX(N-3,I)+    HPCG3370
     1AUX(N-4,I))-.1171875*(AUX(N+6,I)-6.*AUX(N+5,I)-AUX(N+4,I))*H      HPCG3380
     0AUX(N-4,I)=.00390625*(12.*AUX(N-1,I)+135.*AUX(N-2,I)+             HPCG3390
     1108.*AUX(N-3,I)+AUX(N-4,I))-.0234375*(AUX(N+6,I)+18.*AUX(N+5,I)-  HPCG3400
     29.*AUX(N+4,I))*H                                                  HPCG3410
      AUX(N-3,I)=AUX(N-2,I)                                             HPCG3420
  224 AUX(N+4,I)=AUX(N+5,I)                                             HPCG3430
      X=X-H                                                             HPCG3440
      DELT=X-(H+H)                                                      HPCG3450
      CALL FCT(DELT,Y,DERY)                                             HPCG3460
      DO 225 I=1,NDIM                                                   HPCG3470
      AUX(N-2,I)=Y(I)                                                   HPCG3480
      AUX(N+5,I)=DERY(I)                                                HPCG3490
  225 Y(I)=AUX(N-4,I)                                                   HPCG3500
      DELT=DELT-(H+H)                                                   HPCG3510
      CALL FCT(DELT,Y,DERY)                                             HPCG3520
      DO 226 I=1,NDIM                                                   HPCG3530
      DELT=AUX(N+5,I)+AUX(N+4,I)                                        HPCG3540
      DELT=DELT+DELT+DELT                                               HPCG3550
     0AUX(16,I)=8.962963*(AUX(N-1,I)-Y(I))-3.361111*H*(AUX(N+6,I)+DELT  HPCG3560
     1+DERY(I))                                                         HPCG3570
  226 AUX(N+3,I)=DERY(I)                                                HPCG3580
      GOTO 206                                                          HPCG3590
      END                                                               HPCG3600