Google
 

Trailing-Edge - PDP-10 Archives - decuslib10-02 - 43,50145/dhpcg.ssp
There are 2 other files named dhpcg.ssp in the archive. Click here to see a list.
C                                                                       DHCG  10
C     ..................................................................DHCG  20
C                                                                       DHCG  30
C        SUBROUTINE DHPCG                                               DHCG  40
C                                                                       DHCG  50
C        PURPOSE                                                        DHCG  60
C           TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY GENERAL           DHCG  70
C           DIFFERENTIAL EQUATIONS WITH GIVEN INITIAL VALUES.           DHCG  80
C                                                                       DHCG  90
C        USAGE                                                          DHCG 100
C           CALL DHPCG (PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)             DHCG 110
C           PARAMETERS FCT AND OUTP REQUIRE AN EXTERNAL STATEMENT.      DHCG 120
C                                                                       DHCG 130
C        DESCRIPTION OF PARAMETERS                                      DHCG 140
C           PRMT   - DOUBLE PRECISION INPUT AND OUTPUT VECTOR WITH      DHCG 150
C                    DIMENSION GREATER THAN OR EQUAL TO 5, WHICH        DHCG 160
C                    SPECIFIES THE PARAMETERS OF THE INTERVAL AND OF    DHCG 170
C                    ACCURACY AND WHICH SERVES FOR COMMUNICATION BETWEENDHCG 180
C                    OUTPUT SUBROUTINE (FURNISHED BY THE USER) AND      DHCG 190
C                    SUBROUTINE DHPCG. EXCEPT PRMT(5) THE COMPONENTS    DHCG 200
C                    ARE NOT DESTROYED BY SUBROUTINE DHPCG AND THEY ARE DHCG 210
C           PRMT(1)- LOWER BOUND OF THE INTERVAL (INPUT),               DHCG 220
C           PRMT(2)- UPPER BOUND OF THE INTERVAL (INPUT),               DHCG 230
C           PRMT(3)- INITIAL INCREMENT OF THE INDEPENDENT VARIABLE      DHCG 240
C                    (INPUT),                                           DHCG 250
C           PRMT(4)- UPPER ERROR BOUND (INPUT). IF ABSOLUTE ERROR IS    DHCG 260
C                    GREATER THAN PRMT(4), INCREMENT GETS HALVED.       DHCG 270
C                    IF INCREMENT IS LESS THAN PRMT(3) AND ABSOLUTE     DHCG 280
C                    ERROR LESS THAN PRMT(4)/50, INCREMENT GETS DOUBLED.DHCG 290
C                    THE USER MAY CHANGE PRMT(4) BY MEANS OF HIS        DHCG 300
C                    OUTPUT SUBROUTINE.                                 DHCG 310
C           PRMT(5)- NO INPUT PARAMETER. SUBROUTINE DHPCG INITIALIZES   DHCG 320
C                    PRMT(5)=0. IF THE USER WANTS TO TERMINATE          DHCG 330
C                    SUBROUTINE DHPCG AT ANY OUTPUT POINT, HE HAS TO    DHCG 340
C                    CHANGE PRMT(5) TO NON-ZERO BY MEANS OF SUBROUTINE  DHCG 350
C                    OUTP. FURTHER COMPONENTS OF VECTOR PRMT ARE        DHCG 360
C                    FEASIBLE IF ITS DIMENSION IS DEFINED GREATER       DHCG 370
C                    THAN 5. HOWEVER SUBROUTINE DHPCG DOES NOT REQUIRE  DHCG 380
C                    AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL   DHCG 390
C                    FOR HANDING RESULT VALUES TO THE MAIN PROGRAM      DHCG 400
C                    (CALLING DHPCG) WHICH ARE OBTAINED BY SPECIAL      DHCG 410
C                    MANIPULATIONS WITH OUTPUT DATA IN SUBROUTINE OUTP. DHCG 420
C           Y      - DOUBLE PRECISION INPUT VECTOR OF INITIAL VALUES    DHCG 430
C                    (DESTROYED). LATERON Y IS THE RESULTING VECTOR OF  DHCG 440
C                    DEPENDENT VARIABLES COMPUTED AT INTERMEDIATE       DHCG 450
C                    POINTS X.                                          DHCG 460
C           DERY   - DOUBLE PRECISION INPUT VECTOR OF ERROR WEIGHTS     DHCG 470
C                    (DESTROYED). THE SUM OF ITS COMPONENTS MUST BE     DHCG 480
C                    EQUAL TO 1. LATERON DERY IS THE VECTOR OF          DHCG 490
C                    DERIVATIVES, WHICH BELONG TO FUNCTION VALUES Y AT  DHCG 500
C                    INTERMEDIATE POINTS X.                             DHCG 510
C           NDIM   - AN INPUT VALUE, WHICH SPECIFIES THE NUMBER OF      DHCG 520
C                    EQUATIONS IN THE SYSTEM.                           DHCG 530
C           IHLF   - AN OUTPUT VALUE, WHICH SPECIFIES THE NUMBER OF     DHCG 540
C                    BISECTIONS OF THE INITIAL INCREMENT. IF IHLF GETS  DHCG 550
C                    GREATER THAN 10, SUBROUTINE DHPCG RETURNS WITH     DHCG 560
C                    ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM.           DHCG 570
C                    ERROR MESSAGE IHLF=12 OR IHLF=13 APPEARS IN CASE   DHCG 580
C                    PRMT(3)=0 OR IN CASE SIGN(PRMT(3)).NE.SIGN(PRMT(2)-DHCG 590
C                    PRMT(1)) RESPECTIVELY.                             DHCG 600
C           FCT    - THE NAME OF AN EXTERNAL SUBROUTINE USED. IT        DHCG 610
C                    COMPUTES THE RIGHT HAND SIDES DERY OF THE SYSTEM   DHCG 620
C                    TO GIVEN VALUES OF X AND Y. ITS PARAMETER LIST     DHCG 630
C                    MUST BE X,Y,DERY. THE SUBROUTINE SHOULD NOT        DHCG 640
C                    DESTROY X AND Y.                                   DHCG 650
C           OUTP   - THE NAME OF AN EXTERNAL OUTPUT SUBROUTINE USED.    DHCG 660
C                    ITS PARAMETER LIST MUST BE X,Y,DERY,IHLF,NDIM,PRMT.DHCG 670
C                    NONE OF THESE PARAMETERS (EXCEPT, IF NECESSARY,    DHCG 680
C                    PRMT(4),PRMT(5),...) SHOULD BE CHANGED BY          DHCG 690
C                    SUBROUTINE OUTP. IF PRMT(5) IS CHANGED TO NON-ZERO,DHCG 700
C                    SUBROUTINE DHPCG IS TERMINATED.                    DHCG 710
C           AUX    - DOUBLE PRECISION AUXILIARY STORAGE ARRAY WITH 16   DHCG 720
C                    ROWS AND NDIM COLUMNS.                             DHCG 730
C                                                                       DHCG 740
C        REMARKS                                                        DHCG 750
C           THE PROCEDURE TERMINATES AND RETURNS TO CALLING PROGRAM, IF DHCG 760
C           (1) MORE THAN 10 BISECTIONS OF THE INITIAL INCREMENT ARE    DHCG 770
C               NECESSARY TO GET SATISFACTORY ACCURACY (ERROR MESSAGE   DHCG 780
C               IHLF=11),                                               DHCG 790
C           (2) INITIAL INCREMENT IS EQUAL TO 0 OR HAS WRONG SIGN       DHCG 800
C               (ERROR MESSAGES IHLF=12 OR IHLF=13),                    DHCG 810
C           (3) THE WHOLE INTEGRATION INTERVAL IS WORKED THROUGH,       DHCG 820
C           (4) SUBROUTINE OUTP HAS CHANGED PRMT(5) TO NON-ZERO.        DHCG 830
C                                                                       DHCG 840
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  DHCG 850
C           THE EXTERNAL SUBROUTINES FCT(X,Y,DERY) AND                  DHCG 860
C           OUTP(X,Y,DERY,IHLF,NDIM,PRMT) MUST BE FURNISHED BY THE USER.DHCG 870
C                                                                       DHCG 880
C        METHOD                                                         DHCG 890
C           EVALUATION IS DONE BY MEANS OF HAMMINGS MODIFIED PREDICTOR- DHCG 900
C           CORRECTOR METHOD. IT IS A FOURTH ORDER METHOD, USING 4      DHCG 910
C           PRECEEDING POINTS FOR COMPUTATION OF A NEW VECTOR Y OF THE  DHCG 920
C           DEPENDENT VARIABLES.                                        DHCG 930
C           FOURTH ORDER RUNGE-KUTTA METHOD SUGGESTED BY RALSTON IS     DHCG 940
C           USED FOR ADJUSTMENT OF THE INITIAL INCREMENT AND FOR        DHCG 950
C           COMPUTATION OF STARTING VALUES.                             DHCG 960
C           SUBROUTINE DHPCG AUTOMATICALLY ADJUSTS THE INCREMENT DURING DHCG 970
C           THE WHOLE COMPUTATION BY HALVING OR DOUBLING.               DHCG 980
C           TO GET FULL FLEXIBILITY IN OUTPUT, AN OUTPUT SUBROUTINE     DHCG 990
C           MUST BE CODED BY THE USER.                                  DHCG1000
C           FOR REFERENCE, SEE                                          DHCG1010
C           (1)  RALSTON/WILF, MATHEMATICAL METHODS FOR DIGITAL         DHCG1020
C                COMPUTERS, WILEY, NEW YORK/LONDON, 1960, PP.95-109.    DHCG1030
C           (2)  RALSTON, RUNGE-KUTTA METHODS WITH MINIMUM ERROR BOUNDS,DHCG1040
C                MTAC, VOL.16, ISS.80 (1962), PP.431-437.               DHCG1050
C                                                                       DHCG1060
C     ..................................................................DHCG1070
C                                                                       DHCG1080
      SUBROUTINE DHPCG(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)              DHCG1090
C                                                                       DHCG1100
C                                                                       DHCG1110
      DIMENSION PRMT(1),Y(1),DERY(1),AUX(16,1)                          DHCG1120
      DOUBLE PRECISION Y,DERY,AUX,PRMT,X,H,Z,DELT                       DHCG1130
      N=1                                                               DHCG1140
      IHLF=0                                                            DHCG1150
      X=PRMT(1)                                                         DHCG1160
      H=PRMT(3)                                                         DHCG1170
      PRMT(5)=0.D0                                                      DHCG1180
      DO 1 I=1,NDIM                                                     DHCG1190
      AUX(16,I)=0.D0                                                    DHCG1200
      AUX(15,I)=DERY(I)                                                 DHCG1210
    1 AUX(1,I)=Y(I)                                                     DHCG1220
      IF(H*(PRMT(2)-X))3,2,4                                            DHCG1230
C                                                                       DHCG1240
C     ERROR RETURNS                                                     DHCG1250
    2 IHLF=12                                                           DHCG1260
      GOTO 4                                                            DHCG1270
    3 IHLF=13                                                           DHCG1280
C                                                                       DHCG1290
C     COMPUTATION OF DERY FOR STARTING VALUES                           DHCG1300
    4 CALL FCT(X,Y,DERY)                                                DHCG1310
C                                                                       DHCG1320
C     RECORDING OF STARTING VALUES                                      DHCG1330
      CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                DHCG1340
      IF(PRMT(5))6,5,6                                                  DHCG1350
    5 IF(IHLF)7,7,6                                                     DHCG1360
    6 RETURN                                                            DHCG1370
    7 DO 8 I=1,NDIM                                                     DHCG1380
    8 AUX(8,I)=DERY(I)                                                  DHCG1390
C                                                                       DHCG1400
C     COMPUTATION OF AUX(2,I)                                           DHCG1410
      ISW=1                                                             DHCG1420
      GOTO 100                                                          DHCG1430
C                                                                       DHCG1440
    9 X=X+H                                                             DHCG1450
      DO 10 I=1,NDIM                                                    DHCG1460
   10 AUX(2,I)=Y(I)                                                     DHCG1470
C                                                                       DHCG1480
C     INCREMENT H IS TESTED BY MEANS OF BISECTION                       DHCG1490
   11 IHLF=IHLF+1                                                       DHCG1500
      X=X-H                                                             DHCG1510
      DO 12 I=1,NDIM                                                    DHCG1520
   12 AUX(4,I)=AUX(2,I)                                                 DHCG1530
      H=.5D0*H                                                          DHCG1540
      N=1                                                               DHCG1550
      ISW=2                                                             DHCG1560
      GOTO 100                                                          DHCG1570
C                                                                       DHCG1580
   13 X=X+H                                                             DHCG1590
      CALL FCT(X,Y,DERY)                                                DHCG1600
      N=2                                                               DHCG1610
      DO 14 I=1,NDIM                                                    DHCG1620
      AUX(2,I)=Y(I)                                                     DHCG1630
   14 AUX(9,I)=DERY(I)                                                  DHCG1640
      ISW=3                                                             DHCG1650
      GOTO 100                                                          DHCG1660
C                                                                       DHCG1670
C     COMPUTATION OF TEST VALUE DELT                                    DHCG1680
   15 DELT=0.D0                                                         DHCG1690
      DO 16 I=1,NDIM                                                    DHCG1700
   16 DELT=DELT+AUX(15,I)*DABS(Y(I)-AUX(4,I))                           DHCG1710
      DELT=.066666666666666667D0*DELT                                   DHCG1720
      IF(DELT-PRMT(4))19,19,17                                          DHCG1730
   17 IF(IHLF-10)11,18,18                                               DHCG1740
C                                                                       DHCG1750
C     NO SATISFACTORY ACCURACY AFTER 10 BISECTIONS. ERROR MESSAGE.      DHCG1760
   18 IHLF=11                                                           DHCG1770
      X=X+H                                                             DHCG1780
      GOTO 4                                                            DHCG1790
C                                                                       DHCG1800
C     THERE IS SATISFACTORY ACCURACY AFTER LESS THAN 11 BISECTIONS.     DHCG1810
   19 X=X+H                                                             DHCG1820
      CALL FCT(X,Y,DERY)                                                DHCG1830
      DO 20 I=1,NDIM                                                    DHCG1840
      AUX(3,I)=Y(I)                                                     DHCG1850
   20 AUX(10,I)=DERY(I)                                                 DHCG1860
      N=3                                                               DHCG1870
      ISW=4                                                             DHCG1880
      GOTO 100                                                          DHCG1890
C                                                                       DHCG1900
   21 N=1                                                               DHCG1910
      X=X+H                                                             DHCG1920
      CALL FCT(X,Y,DERY)                                                DHCG1930
      X=PRMT(1)                                                         DHCG1940
      DO 22 I=1,NDIM                                                    DHCG1950
      AUX(11,I)=DERY(I)                                                 DHCG1960
   220Y(I)=AUX(1,I)+H*(.375D0*AUX(8,I)+.7916666666666667D0*AUX(9,I)     DHCG1970
     1-.20833333333333333D0*AUX(10,I)+.041666666666666667D0*DERY(I))    DHCG1980
   23 X=X+H                                                             DHCG1990
      N=N+1                                                             DHCG2000
      CALL FCT(X,Y,DERY)                                                DHCG2010
      CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                DHCG2020
      IF(PRMT(5))6,24,6                                                 DHCG2030
   24 IF(N-4)25,200,200                                                 DHCG2040
   25 DO 26 I=1,NDIM                                                    DHCG2050
      AUX(N,I)=Y(I)                                                     DHCG2060
   26 AUX(N+7,I)=DERY(I)                                                DHCG2070
      IF(N-3)27,29,200                                                  DHCG2080
C                                                                       DHCG2090
   27 DO 28 I=1,NDIM                                                    DHCG2100
      DELT=AUX(9,I)+AUX(9,I)                                            DHCG2110
      DELT=DELT+DELT                                                    DHCG2120
   28 Y(I)=AUX(1,I)+.33333333333333333D0*H*(AUX(8,I)+DELT+AUX(10,I))    DHCG2130
      GOTO 23                                                           DHCG2140
C                                                                       DHCG2150
   29 DO 30 I=1,NDIM                                                    DHCG2160
      DELT=AUX(9,I)+AUX(10,I)                                           DHCG2170
      DELT=DELT+DELT+DELT                                               DHCG2180
   30 Y(I)=AUX(1,I)+.375D0*H*(AUX(8,I)+DELT+AUX(11,I))                  DHCG2190
      GOTO 23                                                           DHCG2200
C                                                                       DHCG2210
C     THE FOLLOWING PART OF SUBROUTINE DHPCG COMPUTES BY MEANS OF       DHCG2220
C     RUNGE-KUTTA METHOD STARTING VALUES FOR THE NOT SELF-STARTING      DHCG2230
C     PREDICTOR-CORRECTOR METHOD.                                       DHCG2240
  100 DO 101 I=1,NDIM                                                   DHCG2250
      Z=H*AUX(N+7,I)                                                    DHCG2260
      AUX(5,I)=Z                                                        DHCG2270
  101 Y(I)=AUX(N,I)+.4D0*Z                                              DHCG2280
C     Z IS AN AUXILIARY STORAGE LOCATION                                DHCG2290
C                                                                       DHCG2300
      Z=X+.4D0*H                                                        DHCG2310
      CALL FCT(Z,Y,DERY)                                                DHCG2320
      DO 102 I=1,NDIM                                                   DHCG2330
      Z=H*DERY(I)                                                       DHCG2340
      AUX(6,I)=Z                                                        DHCG2350
  102 Y(I)=AUX(N,I)+.29697760924775360D0*AUX(5,I)+.15875964497103583D0*ZDHCG2360
C                                                                       DHCG2370
      Z=X+.45573725421878943D0*H                                        DHCG2380
      CALL FCT(Z,Y,DERY)                                                DHCG2390
      DO 103 I=1,NDIM                                                   DHCG2400
      Z=H*DERY(I)                                                       DHCG2410
      AUX(7,I)=Z                                                        DHCG2420
  103 Y(I)=AUX(N,I)+.21810038822592047D0*AUX(5,I)-3.0509651486929308D0* DHCG2430
     1AUX(6,I)+3.8328647604670103D0*Z                                   DHCG2440
C                                                                       DHCG2450
      Z=X+H                                                             DHCG2460
      CALL FCT(Z,Y,DERY)                                                DHCG2470
      DO 104 I=1,NDIM                                                   DHCG2480
  1040Y(I)=AUX(N,I)+.17476028226269037D0*AUX(5,I)-.55148066287873294D0* DHCG2490
     1AUX(6,I)+1.2055355993965235D0*AUX(7,I)+.17118478121951903D0*      DHCG2500
     2H*DERY(I)                                                         DHCG2510
      GOTO(9,13,15,21),ISW                                              DHCG2520
C                                                                       DHCG2530
C     POSSIBLE BREAK-POINT FOR LINKAGE                                  DHCG2540
C                                                                       DHCG2550
C     STARTING VALUES ARE COMPUTED.                                     DHCG2560
C     NOW START HAMMINGS MODIFIED PREDICTOR-CORRECTOR METHOD.           DHCG2570
  200 ISTEP=3                                                           DHCG2580
  201 IF(N-8)204,202,204                                                DHCG2590
C                                                                       DHCG2600
C     N=8 CAUSES THE ROWS OF AUX TO CHANGE THEIR STORAGE LOCATIONS      DHCG2610
  202 DO 203 N=2,7                                                      DHCG2620
      DO 203 I=1,NDIM                                                   DHCG2630
      AUX(N-1,I)=AUX(N,I)                                               DHCG2640
  203 AUX(N+6,I)=AUX(N+7,I)                                             DHCG2650
      N=7                                                               DHCG2660
C                                                                       DHCG2670
C     N LESS THAN 8 CAUSES N+1 TO GET N                                 DHCG2680
  204 N=N+1                                                             DHCG2690
C                                                                       DHCG2700
C     COMPUTATION OF NEXT VECTOR Y                                      DHCG2710
      DO 205 I=1,NDIM                                                   DHCG2720
      AUX(N-1,I)=Y(I)                                                   DHCG2730
  205 AUX(N+6,I)=DERY(I)                                                DHCG2740
      X=X+H                                                             DHCG2750
  206 ISTEP=ISTEP+1                                                     DHCG2760
      DO 207 I=1,NDIM                                                   DHCG2770
     0DELT=AUX(N-4,I)+1.3333333333333333D0*H*(AUX(N+6,I)+AUX(N+6,I)-    DHCG2780
     1AUX(N+5,I)+AUX(N+4,I)+AUX(N+4,I))                                 DHCG2790
      Y(I)=DELT-.9256198347107438D0*AUX(16,I)                           DHCG2800
  207 AUX(16,I)=DELT                                                    DHCG2810
C     PREDICTOR IS NOW GENERATED IN ROW 16 OF AUX, MODIFIED PREDICTOR   DHCG2820
C     IS GENERATED IN Y. DELT MEANS AN AUXILIARY STORAGE.               DHCG2830
C                                                                       DHCG2840
      CALL FCT(X,Y,DERY)                                                DHCG2850
C     DERIVATIVE OF MODIFIED PREDICTOR IS GENERATED IN DERY             DHCG2860
C                                                                       DHCG2870
      DO 208 I=1,NDIM                                                   DHCG2880
     0DELT=.125D0*(9.D0*AUX(N-1,I)-AUX(N-3,I)+3.D0*H*(DERY(I)+AUX(N+6,I)DHCG2890
     1+AUX(N+6,I)-AUX(N+5,I)))                                          DHCG2900
      AUX(16,I)=AUX(16,I)-DELT                                          DHCG2910
  208 Y(I)=DELT+.07438016528925620D0*AUX(16,I)                          DHCG2920
C                                                                       DHCG2930
C     TEST WHETHER H MUST BE HALVED OR DOUBLED                          DHCG2940
      DELT=0.D0                                                         DHCG2950
      DO 209 I=1,NDIM                                                   DHCG2960
  209 DELT=DELT+AUX(15,I)*DABS(AUX(16,I))                               DHCG2970
      IF(DELT-PRMT(4))210,222,222                                       DHCG2980
C                                                                       DHCG2990
C     H MUST NOT BE HALVED. THAT MEANS Y(I) ARE GOOD.                   DHCG3000
  210 CALL FCT(X,Y,DERY)                                                DHCG3010
      CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                DHCG3020
      IF(PRMT(5))212,211,212                                            DHCG3030
  211 IF(IHLF-11)213,212,212                                            DHCG3040
  212 RETURN                                                            DHCG3050
  213 IF(H*(X-PRMT(2)))214,212,212                                      DHCG3060
  214 IF(DABS(X-PRMT(2))-.1D0*DABS(H))212,215,215                       DHCG3070
  215 IF(DELT-.02D0*PRMT(4))216,216,201                                 DHCG3080
C                                                                       DHCG3090
C                                                                       DHCG3100
C     H COULD BE DOUBLED IF ALL NECESSARY PRECEEDING VALUES ARE         DHCG3110
C     AVAILABLE                                                         DHCG3120
  216 IF(IHLF)201,201,217                                               DHCG3130
  217 IF(N-7)201,218,218                                                DHCG3140
  218 IF(ISTEP-4)201,219,219                                            DHCG3150
  219 IMOD=ISTEP/2                                                      DHCG3160
      IF(ISTEP-IMOD-IMOD)201,220,201                                    DHCG3170
  220 H=H+H                                                             DHCG3180
      IHLF=IHLF-1                                                       DHCG3190
      ISTEP=0                                                           DHCG3200
      DO 221 I=1,NDIM                                                   DHCG3210
      AUX(N-1,I)=AUX(N-2,I)                                             DHCG3220
      AUX(N-2,I)=AUX(N-4,I)                                             DHCG3230
      AUX(N-3,I)=AUX(N-6,I)                                             DHCG3240
      AUX(N+6,I)=AUX(N+5,I)                                             DHCG3250
      AUX(N+5,I)=AUX(N+3,I)                                             DHCG3260
      AUX(N+4,I)=AUX(N+1,I)                                             DHCG3270
      DELT=AUX(N+6,I)+AUX(N+5,I)                                        DHCG3280
      DELT=DELT+DELT+DELT                                               DHCG3290
  2210AUX(16,I)=8.962962962962963D0*(Y(I)-AUX(N-3,I))                   DHCG3300
     1-3.3611111111111111D0*H*(DERY(I)+DELT+AUX(N+4,I))                 DHCG3310
      GOTO 201                                                          DHCG3320
C                                                                       DHCG3330
C                                                                       DHCG3340
C     H MUST BE HALVED                                                  DHCG3350
  222 IHLF=IHLF+1                                                       DHCG3360
      IF(IHLF-10)223,223,210                                            DHCG3370
  223 H=.5D0*H                                                          DHCG3380
      ISTEP=0                                                           DHCG3390
      DO 224 I=1,NDIM                                                   DHCG3400
     0Y(I)=.390625D-2*(8.D1*AUX(N-1,I)+135.D0*AUX(N-2,I)+4.D1*AUX(N-3,I)DHCG3410
     1+AUX(N-4,I))-.1171875D0*(AUX(N+6,I)-6.D0*AUX(N+5,I)-AUX(N+4,I))*H DHCG3420
     0AUX(N-4,I)=.390625D-2*(12.D0*AUX(N-1,I)+135.D0*AUX(N-2,I)+        DHCG3430
     1108.D0*AUX(N-3,I)+AUX(N-4,I))-.0234375D0*(AUX(N+6,I)+             DHCG3440
     218.D0*AUX(N+5,I)-9.D0*AUX(N+4,I))*H                               DHCG3450
      AUX(N-3,I)=AUX(N-2,I)                                             DHCG3460
  224 AUX(N+4,I)=AUX(N+5,I)                                             DHCG3470
      X=X-H                                                             DHCG3480
      DELT=X-(H+H)                                                      DHCG3490
      CALL FCT(DELT,Y,DERY)                                             DHCG3500
      DO 225 I=1,NDIM                                                   DHCG3510
      AUX(N-2,I)=Y(I)                                                   DHCG3520
      AUX(N+5,I)=DERY(I)                                                DHCG3530
  225 Y(I)=AUX(N-4,I)                                                   DHCG3540
      DELT=DELT-(H+H)                                                   DHCG3550
      CALL FCT(DELT,Y,DERY)                                             DHCG3560
      DO 226 I=1,NDIM                                                   DHCG3570
      DELT=AUX(N+5,I)+AUX(N+4,I)                                        DHCG3580
      DELT=DELT+DELT+DELT                                               DHCG3590
     0AUX(16,I)=8.962962962962963D0*(AUX(N-1,I)-Y(I))                   DHCG3600
     1-3.3611111111111111D0*H*(AUX(N+6,I)+DELT+DERY(I))                 DHCG3610
  226 AUX(N+3,I)=DERY(I)                                                DHCG3620
      GOTO 206                                                          DHCG3630
      END                                                               DHCG3640