      SUBROUTINE MUSN(FDIF,X0T,G,N,A,B,ER,TI,NTI,NRTI,AMP,ITLIM,X,Q,
     1                U,NU,D,PHI,KP,W,LW,IW,LIW,WGR,LWG,IERROR)
C     ---------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ER(5),TI(NTI),X(N,NTI),Q(N,N,NTI),U(NU,NTI),D(N,NTI),
     1          PHI(NU,NTI),W(LW),ER1(5),WGR(LWG)
      INTEGER IW(LIW)
      CHARACTER DOOR
      EXTERNAL FDIF,X0T,G
      DIMENSION ALINC(4)
      LOGICAL DIAGNO
C
      DIAGNO = IERROR.EQ.1
      IERROR = 0
C     CHECK INPUT PARAMETERS
C     INPUT ERROR 105
      IF ((N.LT.1).OR.(NRTI.LT.0).OR.(NTI.LT.3).OR.(NU.LT.N*(N+1)/2)
     1 .OR.(A.EQ.B)) GOTO 5000
C     INPUT ERROR 101
      IF ((ER(1).LT.0.D0).OR.(ER(2).LT.0.D0).OR.(ER(3).LT.0.D0))
     1 GOTO 5010
C     INPUT ERROR 106
      K = N * NTI
      L = N * N
      IF ((LW.LT.(7*N+3*K+4*L)).OR.(LIW.LT.3*N+NTI)) GOTO 5020
C     SETTING POINTERS FOR THE WORK ARRAYS W AND IW
      I1 = 7 * N + 1
      I2 = I1 + K
      I3 = I2 + K
      I4 = I3 + K
      I5 = I4 + L
      I6 = I5 + L
      I7 = I6 + L
      J1 = N + 1
      J2 = J1 + N
      J3 = J2 + N
      TOL = ER(1)
      ER1(3) = ER(3)
      ER1(1) = 1.1D-12 + 3.D0 * ER(3)
      TOL1 = DMAX1(ER(1),DMIN1(ER(2),1.D-2))
      ER1(2) = TOL1
      EPS = DSQRT(ER(3))
      ALINCR = DSQRT(ER(1)/ER(3)) / 4.D0
      IF ((AMP.LE.1.D0).OR.(AMP.GT.ALINCR)) THEN
        ALINC(1) = ALINCR
      ELSE
        ALINC(1) = AMP
      ENDIF
      ALINC(2) = ALINC(1) * ALINC(1)
      ALINC(3) = ALINC(1) * EPS
      ALINC(4) = ALINC(3) * ALINC(1)
      IF (DIAGNO) WRITE(*,200) TOL1
      ITER = 0
      ITTOT = 0
      IHL = 0
      ALAMIN = 0.007D0
 1000 CALL DCSHPO(FDIF,X0T,N,A,B,ER1,TI,NTI,NRTI,ALINC,X,W(I1),Q,U,NU,
     1            W(1),W(I2),W(I4),W(I5),IW(1),IW(J1),WGR,LWG,IERROR)
      NEG = IW(J1)
      IF ((IERROR.NE.0).AND.(IERROR.NE.213)) GOTO 6000
      IF (IERROR.NE.0) THEN
        WRITE(*,110) IERROR
        CALL ERRHAN(IERROR,ER,NEG)
      ENDIF
      NOSH = NRTI
      JAC = 1
      NG = 0
      CALL DCSAOJ(FDIF,N,ER1,TI,NTI,NRTI,X,JAC,NG,IW(1),W(I1),Q,U,NU,
     1            W(1),W(I4),W(I5),IW(J1),WGR,LWG,IERROR)
      NEG = IW(J1)
      IF ((IERROR.NE.0).AND.(IERROR.NE.213)) GOTO 6000
      IF (DIAGNO) THEN
        WRITE(*,210)
        DO 1050 I = 1 , NRTI , 5
          II = MIN0(NRTI,I+4)
          WRITE(*,220) (TI(J),J=I,II)
 1050   CONTINUE
      ENDIF
      JN = I1 + N * NRTI -1
 1100 CALL DJINGX(G,N,NRTI,X,W(I1),NTI,Q,U,NU,ER,IW(1),JAC,D,KP,PHI,
     1            W(1),W(I4),W(I5),W(I6),W(I7),IW(J1),IW(J2),IERROR)
      IF ((IERROR.NE.0).AND.(IERROR.NE.240)) GOTO 6000
      IF (IERROR.NE.0) THEN
        WRITE(*,110) IERROR
        CALL ERRHAN(IERROR,ER,NEG)
      ENDIF
C     XI IN W(I1)
      NM = N * NRTI
      CALL DINPRO(W(I1),1,W(I1),1,NM,RKSI)
      RKSI = DSQRT(RKSI)
      IF (DIAGNO) THEN
        WRITE(*,230) RKSI
        WRITE(*,240) ER(4),ER(5)
      ENDIF
      IF (RKSI.LT.TOL1) GOTO 2500
      ALAM = 1.D0
      ALM1 = 1.D0
C     COMPUTE NEXT ITERATION
 1400 DO 1500 I = 1 , NRTI
        L1 = (I-1) * N + I3 - 1
        L2 = (I-1) * N + I1 - 1
        DO 1500 J = 1 , N
          W(L1+J) = X(J,I) + ALAM * W(L2+J)
 1500 CONTINUE
C     SI+ STORED IN W(I3)
      IF (DIAGNO) THEN
        WRITE(*,250) ALAM
      ENDIF
      NG = 0
      CALL DCSAOJ(FDIF,N,ER1,TI,NTI,NRTI,W(I3),JAC,NG,IW(1),W(I2),Q,U,
     1            NU,W(1),W(I4),W(I5),IW(J1),WGR,LWG,IERROR)
      NEG = IW(J1)
      IF ((IERROR.NE.0).AND.(IERROR.NE.213)) GOTO 6000
      IF (IERROR.NE.0) THEN
        WRITE(*,110) IERROR
        CALL ERRHAN(IERROR,ER,NEG)
      ENDIF
      IF (JAC.NE.0) CALL DCHINC(W(I3),W(I2),N,NTI,U,NU,NRTI,ALINC(2),
     1                          IW(1),IW(J3))
      CALL DJINGX(G,N,NRTI,W(I3),W(I2),NTI,Q,U,NU,ER,IW(1),JAC,D,KP,PHI,
     1            W(1),W(I4),W(I5),W(I6),W(I7),IW(J1),IW(J2),IERROR)
C     XI+ STORED IN W(I2)
      IF ((IERROR.EQ.260).AND.(JAC.NE.0)) THEN
        ALAM = ALAM / 2.D0
        IF (DIAGNO) THEN
          WRITE(*,260) ALAM
        ENDIF
        IF (ALAM.LT.ALAMIN) GOTO 5100
        CALL DCSAOJ(FDIF,N,ER1,TI,NTI,NRTI,X,JAC,NG,IW(1),W(I1),Q,U,NU,
     1              W(1),W(I4),W(I5),IW(J1),WGR,LWG,IERROR)
        NEG = IW(J1)
        CALL DJINGX(G,N,NRTI,X,W(I1),NTI,Q,U,NU,ER,IW(1),JAC,D,KP,PHI,
     1              W(1),W(I4),W(I5),W(I6),W(I7),IW(J1),IW(J2),IERROR)
        JAC = 0
        IW(J3) = 0
        GOTO 1400
      ENDIF
      IF ((IERROR.NE.0).AND.(IERROR.NE.240)) GOTO 6000
      IF (IERROR.NE.0) THEN
        WRITE(*,110) IERROR
        CALL ERRHAN(IERROR,ER,NEG)
      ENDIF
      CALL DINPRO(W(I2),1,W(I2),1,NM,RKSI0)
      RKSI0 = DSQRT(RKSI0)
      IF (DIAGNO) THEN
        WRITE(*,230) RKSI0
        WRITE(*,240) ER(4),ER(5)
      ENDIF
      IF (RKSI0.GT.RKSI) THEN
        ALAM = ALAM / 2.D0
        IHL = 1
        IF (ALAM.LT.ALAMIN) GOTO 3000
        IF (JAC.EQ.1) THEN
        CALL DCSAOJ(FDIF,N,ER1,TI,NTI,NRTI,X,JAC,NG,IW(1),W(I1),Q,U,
     1              NU,W(1),W(I4),W(I5),IW(J1),WGR,LWG,IERROR)
        NEG = IW(J1)
        CALL DJINGX(G,N,NRTI,X,W(I1),NTI,Q,U,NU,ER,IW(1),JAC,D,KP,PHI,
     1              W(1),W(I4),W(I5),W(I6),W(I7),IW(J1),IW(J2),IERROR)
        ENDIF
        JAC = 0
        IW(J3) = 0
        GOTO 1400
      ELSE
        JAC1 = JAC
        IF ((RKSI0.LT.0.6D0*RKSI).AND.(ALAM.EQ.1.D0)) JAC = 1
        IF (JAC1.NE.JAC) THEN
          DO 1650 I = J3 , J3+NRTI-1
             IW(I) = 0
 1650     CONTINUE
          IW(J3) = -1
        ENDIF
        RKSI = RKSI0
      ENDIF
C     LAMBDA IS ACCEPTED
C     COMPUTING NEXT PREDICTED LAMBDA
      IF ((ALAM.EQ.ALM1).AND.(IHL.EQ.0)) THEN
        ALAM = DMIN1(2.D0 * ALAM,1.D0)
      ELSE
        IHL = 0
        ALM1 = ALAM
      ENDIF
      DO 1700 I = 1 , NRTI
        L1 = (I-1) * N + I3 - 1
        DO 1700 J = 1 , N
          X(J,I) = W(L1+J)
 1700 CONTINUE
C     X := SI+ STORED IN X
      DO 1800 I = 1 , NM
        W(I1-1+I) = W(I2-1+I)
 1800 CONTINUE
C     XI := XI+  STORED IN W(I1)
      ITER = ITER + 1
      ITTOT = ITTOT + 1
      IF (DIAGNO) THEN
        IF (ALM1.EQ.1.D0) THEN
          WRITE(*,280) ITTOT
        ELSE
          WRITE(*,290) ITTOT,ALM1
        ENDIF
      ENDIF
      IF (RKSI.LT.TOL1) GOTO 2500
      IF (ITER.GT.ITLIM) GOTO 5200
      IF ((JAC.EQ.0).OR.(IW(J3).EQ.0)) GOTO 1400
 2400 CALL DNEWPO(FDIF,N,TOL1,ER(3),TI,NTI,NRTI,ALINC,X,W(I1),IW(1),Q,
     1            U,NU,IW(J3),WGR,LWG,W(1),W(I4),W(I5),IERROR)
      NEG = IW(J3)
      IF (IERROR.NE.0) GOTO 6000
      IF (DIAGNO.AND.NRTI.NE.NOSH) THEN
        WRITE(*,300) NRTI-NOSH
        NOSH = NRTI
      ENDIF
      GOTO 1100
 2500 IF (TOL1.LE.TOL) THEN
       WRITE(*,*) 'RETURN FROM MUSN : NO ITERATIONS =',ITTOT
       RETURN
      ENDIF
      TOL1 = DMAX1(TOL1*TOL1,TOL)
      IF (DABS(TOL-TOL1).LT.100.D0*ER(3)) TOL1 = TOL
      ER1(1) = 1.1D-12 + 3.D0 * ER(3)
      ER1(2) = TOL1
      ITER = 0
      IF (DIAGNO) THEN
        WRITE(*,310) TOL1
      ENDIF
      JAC = 0
      NG = 1
      CALL DCSAOJ(FDIF,N,ER1,TI,NTI,NRTI,X,JAC,NG,IW(1),W(I1),Q,U,NU,
     1            W(1),W(I4),W(I5),IW(J1),WGR,LWG,IERROR)
      NEG = IW(J1)
      IF ((IERROR.NE.0).AND.(IERROR.NE.213)) GOTO 6000
      IW(J3) = -1
      DO 2600 I = J3+1 , J3+NRTI-1
        IW(I) = 0
 2600 CONTINUE
      JAC = 1
      GOTO 2400
 3000 TOL1 = 1.D-1 * TOL1
      IF (DABS(TOL-TOL1).LT.100.D0*ER(3)) TOL1 = TOL
      ER1(1) = 1.1D-12 + 3.D0 * ER(3)
      ER1(2) = TOL1
      IF (TOL1.LT.TOL) GOTO 5100
      IF (DIAGNO) THEN
        WRITE(*,320) TOL1
      ENDIF
      NRTI = 1
      GOTO 1000
 5000 IERROR = 105
      GOTO 6000
 5010 IERROR = 101
      GOTO 6000
 5020 IERROR = 106
      GOTO 6000
 5100 IERROR = 230
      GOTO 6000
 5200 IERROR = 231
 6000 WRITE(*,100) IERROR
      CALL ERRHAN(IERROR,ER,NEG)
      RETURN
  100 FORMAT(' TERMINAL ERROR IN MUSN : IERROR = ',I4,/)
  110 FORMAT(' WARNING ERROR IN MUSN : IERROR = ',I4,/)
  200 FORMAT(' ',4('*****'),' START TOLERANCE : ',8X,1P,D12.5)
  210 FORMAT(5X,'SHOOTING POINTS :')
  220 FORMAT(5X,5(D12.5,3X))
  230 FORMAT(5X,'NORMALIZED RESIDUE : ',22X,1P,D16.9)
  240 FORMAT(5X,'COND.NUMBER = ',1P,D12.5,2X,'AMPLI.FACTOR = ',D12.5)
  250 FORMAT(5X,'PREDICTED DAMPING FACTOR : ',16X,F10.7)
  260 FORMAT(' IERROR = 260 , NEWTON UPDATE CAUSES ILL-CONDITIONING',/,
     1 ' WITH RESPECT TO THE BC; WILL TRY A DAMPING FACTOR=',F10.7)
  280 FORMAT(' ',10('******'),/,' ',4('****'),2X,I4,
     1 ' ITERATION COMPLETED',2X,4('****'),/,' ',4('*****'),2X,
     2 'FULL NEWTON STEP',2X,4('*****'),/,' ',10('******'),/)
  290 FORMAT(' ',10('******'),/,' ',4('****'),2X,I4,
     1 ' ITERATION COMPLETED',2X,4('****'),/,' ***** ',
     2 'DAMPED NEWTON STEP ; DAMPING FACTOR = ',F10.7,' *****',/,' ',
     3 10('******'),/)
  300 FORMAT(' ',I3,' NEW SHOOTING POINTS INSERTED.',/)
  310 FORMAT(' ',4('*****'),' NEW TOLERANCE : ',10X,1P,D12.5)
  320 FORMAT(' ***** NEWTON FAILED WITH THIS TOLERANCE. ***** ',/,
     1 ' ***** WILL TRY NEW TOLERANCE : ',10X,'TOL = ',1P,D12.5)
C     END OF MUSN
      END
      SUBROUTINE DCHINC(X,S,N,NTI,U,NU,NRTI,ALINCR,KKK,IPG)
C     -----------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N,NTI),S(N,NTI),U(NU,NTI)
      INTEGER KKK(N),IPG(NTI)
C
      IPG(1) = 0
      DO 1500 K = 2 , NRTI
        KM1 = K - 1
        E1 = 0.D0
        E2 = 0.D0
        DO 1000 J = 1 , N
          E3 = DABS(X(J,KM1))
          IF (E3.GT.E1) E1 = E3
          E3 = DABS(S(J,K))
          IF (E3.GT.E2) E2 = E3
 1000   CONTINUE
        E3 = E2 / E1
        IF (E3.GT.ALINCR) THEN
          IPG(1) = 1
          IPG(K) = 1
        ELSE
          IPG(K) = 0
        ENDIF
 1500 CONTINUE
      DO 2500 K = 2 , NRTI
        E1 = 0.D0
        DO 2000 J = 1 , N
          E2 = U(KKK(J),K)
          IF (E2.GT.E1) E1 = E2
 2000   CONTINUE
        IF (E1.GT.ALINCR) THEN
          IPG(1) = 1
          IPG(K) = 1
        ENDIF
 2500 CONTINUE
      RETURN
C     END OF DCHINC
      END
      SUBROUTINE DCSHPO(FDIF,X0T,N,A,B,ER,TI,NTI,NRTI,ALINC,X,S,Q,U,NU,
     1                  W,WTI,WF,WF0,KKK,IP,WGR,LWGR,IERROR)
C     ----------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ER(5),TI(NTI),ALINC(4),X(N,NTI),S(N,NTI),Q(N,N,NTI),
     1          U(NU,NTI),W(N,7),WTI(NTI),WF(N,N),WF0(N,N),WGR(LWGR)
      INTEGER KKK(N),IP(N)
      EXTERNAL FDIF,X0T
C
      IERROR = 0
      EPS = DSQRT(ER(3))
      TI(1) = A
      WTI(1) = A
      NOTI = 2
      IF (NRTI.EQ.0) TI(2) = B
      IF (NRTI.EQ.1) THEN
        IF (A.LE.B) THEN
          DO 1100 I = 2 , NTI
            IF (TI(I-1).GE.TI(I)) GOTO 5200
            IF (TI(I).GE.B) GOTO 1200
 1100     CONTINUE
 1200     NOTI=I
        ELSE
          DO 1300 I = 2 , NTI
            IF (TI(I-1).LE.TI(I)) GOTO 5200
            IF (TI(I).LE.B) GOTO 1400
 1300     CONTINUE
 1400     NOTI=I
        ENDIF
        IF (TI(NOTI).NE.B) GOTO 5300
      ENDIF
      IF (NRTI.GT.1) THEN
        NOTI = NRTI + 1
        ST = (B-A) / NRTI
        TI(NRTI+1) = B
        DO 1500 I = 2 , NRTI
          TI(I) = TI(I-1) + ST
 1500   CONTINUE
      ENDIF
      IF (LWGR.LT.NOTI) GOTO 5550
C     DETERMINATION OF THE SHOOTING POINTS. GIVEN OUTPUTPOINTS WILL
C     BE SHOOTING POINTS. IF THE INCREMENT BETWEEN TWO OUTPUTPOINTS
C     IS GREATER THAN ALINC(1), A NEW SHOOTING POINT IS INSERTED.
C     DURING THE DETERMINATION THE SHOOTING POINTS WILL BE STORED IN
C     THE ARRAY WTI.
      DO 1600 I = 1 , N
        KKK(I) = I*(I+1) / 2
      DO 1600 J = 1 , N
        IF (I.EQ.J) THEN
          Q(I,I,1) = 1.D0
        ELSE
          Q(I,J,1) = 0.D0
        ENDIF
        Q(I,J,NTI) = Q(I,J,1)
 1600 CONTINUE
      WTI(1) = TI(1)
      JTI = 2
      NHI = 5
      NRTI = 2
      MWGR = 0
      DO 2600 K = 2 , NTI
        KM1 = K - 1
        T1 = WTI(KM1)
        T2 = TI(JTI)
 1650   CALL X0T(T1,X(1,KM1))
        IF (T1.EQ.B) GOTO 2700
        X1 = 0.D0
        DO 1700 I = 1 , N
          ST = DABS(X(I,KM1))
          S(I,K) = X(I,KM1)
          IF (ST.GT.X1) X1 = ST
          DO 1700 J = 1 , N
            WF(J,I) = X(J,KM1) + EPS * Q(J,I,KM1)
 1700   CONTINUE
        IF (X1.EQ.0.D0) X1 = 1.D0
        IWGR = 0
        IER = 1
        IF (K.GT.2) IER=-1
 1900   CALL DRKFMS(FDIF,N,S(1,K),T1,T2,ER(1),ER(2),ER(3),WF,N,HI,NHI,
     1              W,IER)
        IF (IER.GT.3) GOTO 5400
C     CHECK WHETHER T2 SHOULD BE A SHOOTING POINT
        X2 = 0.D0
        X3 = 0.D0
        DO 2000 I = 1 , N
          ST = DABS(S(I,K))
          IF (ST.GT.X2) X2 = ST
          DO 2000 J = 1 , N
            ST = DABS(WF(J,I)-S(J,K))
            IF (ST.GT.X3) X3 = ST
 2000   CONTINUE
        IF ((X2/X1.LT.ALINC(1)).AND.(X3.LT.ALINC(3))) THEN
          IWGR = IWGR + 1
          MWGR = MWGR + 1
          IF (MWGR.GT.LWGR) GOTO 5600
          T1 = T2
          WGR(MWGR) = T2
          IER = -1
          IF (T2.EQ.TI(JTI)) GOTO 2040
          T2 = TI(JTI)
          GOTO 1900
        ENDIF
        IF ((X2/X1.LE.ALINC(2)).AND.(X3.LE.ALINC(4))) THEN
          MWGR = MWGR + 1
          IF (MWGR.GT.LWGR) GOTO 5600
          IWGR = IWGR + 1
          WGR(MWGR) = T2
          IF (T2.EQ.TI(JTI)) GOTO 2040
          WTI(K) = T2
          GOTO 2050
        ENDIF
C     INCREMENT GREATER THAN ALINC(2)
        IF (IWGR.EQ.0) THEN
          T1 = WTI(KM1)
          T2 = (T2 + T1) / 2.D0
          GOTO 1650
        ELSE
          WTI(K) = WGR(MWGR)
          GOTO 2050
        ENDIF
C     A SHOOTING POINT IS DETERMINED.
C     COMPUTATION OF JK, QK AND UK
 2040   JTI = JTI + 1
        WTI(K) = T2
 2050   DO 2100 I = 1 , N
        DO 2100 J = 1 , N
          WF(J,I) = (WF(J,I) - S(J,K)) / EPS
          IF (K.EQ.2) WF0(J,I) = WF(J,I)
 2100   CONTINUE
        CALL DQUDEC(WF,N,N,N,N,W(1,1),W(1,2),W(1,3))
        IF (K.EQ.2) THEN
          ISORT = 0
          DO 2400 ITEL = 1 , N
            CALL DSORTD(W(1,1),N,N,KP,IP,ISORT,IER)
            IF (IER.EQ.0) GOTO 2500
            DO 2200 I = 1 , N
              IK = IP(I)
            DO 2200 J = 1 , N
              WF(J,I) = WF0(J,IK)
              Q(J,I,1) = Q(J,IK,NTI)
 2200       CONTINUE
            DO 2300 I = 1 , N
            DO 2300 J = 1 , N
              WF0(I,J) = WF(I,J)
              Q(I,J,NTI) = Q(I,J,1)
 2300       CONTINUE
            CALL DQUDEC(WF,N,N,N,N,W(1,1),W(1,2),W(1,3))
 2400     CONTINUE
        ENDIF
 2500   CALL DFQUS(WF,N,N,N,W(1,1),W(1,2),K,KKK,U,NU,NTI,Q,WF0,W(1,3))
        NRTI = K
 2600 CONTINUE
 2700 DO 2800 I = 1 , NRTI
        TI(I) = WTI(I)
 2800 CONTINUE
      IF (TI(NRTI).NE.B) GOTO 5500
      RETURN
 5000 IERROR = 100
      RETURN
 5100 IERROR = 101
      RETURN
 5200 IERROR = 120
      RETURN
 5300 IERROR = 121
      RETURN
 5400 IERROR = 210 + IER
      RETURN
 5500 IERROR = 122
      RETURN
 5550 IERROR = 123
      IP(1) = NOTI
      RETURN
 5600 IERROR = 219
      IP(1) = (B - A) / (T1 - A) * (LWGR+1) + 1
      RETURN
C     END OF DCSHPO
      END
      SUBROUTINE DCSAOJ(FDIF,N,ER,TI,NTI,NRTI,X,JA,NG,KKK,S,Q,U,NU,W,
     1                  WF,WF0,IP,WGR,LWGR,IERROR)
C     ----------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ER(5),TI(NTI),X(N,NTI),S(N,NTI),Q(N,N,NTI),U(NU,NTI),
     1          W(N,7),WF(N,N),WF0(N,N),WGR(LWGR)
      INTEGER KKK(N),IP(N)
      DIMENSION HS(5)
      EXTERNAL FDIF
C
C     INPUT ERROR 122
      IF (NTI-NRTI.LT.1) GOTO 5000
      IERROR = 0
      NHI = 0
      EPS = DSQRT(ER(3))
      IF (JA.EQ.1) THEN
        DO 1100 I = 1 , N
        DO 1100 J = 1 , N
          Q(I,J,NTI) = Q(I,J,1)
 1100   CONTINUE
      ENDIF
      MWGR = 1
      IF (NG.EQ.1) MWGR = 0
      DO 1900 K = 2 , NRTI
        KM1 = K - 1
        T1 = TI(KM1)
        T2 = TI(K)
        DO 1200 I = 1 , N
          S(I,K) = X(I,KM1)
 1200   CONTINUE
        IF (JA.EQ.1) THEN
          DO 1300 I = 1 , N
          DO 1300 J = 1 , N
            WF(J,I) = X(J,KM1) + EPS * Q(J,I,KM1)
 1300     CONTINUE
          NF = N
        ELSE
          NF = 0
        ENDIF
        IF (NG.EQ.1) THEN
         NHI = 5
         IER = -1
         IF (K.EQ.2) IER = 1
 1320    CALL DRKFMS(FDIF,N,S(1,K),T1,T2,ER(1),ER(2),ER(3),WF,NF,HI,NHI,
     1               W,IER)
         IF (IER.GT.3) GOTO 5100
         MWGR = MWGR + 1
         IF (MWGR.GT.LWGR) GOTO 5200
         WGR(MWGR) = T2
         IF (T2.NE.TI(K)) THEN
           T1 = T2
           T2 = TI(K)
           IER = -1
           GOTO 1320
         ENDIF
        ELSE
          T0 = T1
 1340     ST = (WGR(MWGR) - T0) / 5.D0
          DO 1360 I = 1 , 5
            HS(I) = ST
 1360     CONTINUE
          CALL DRKFGS(FDIF,N,S(1,K),T1,HS,5,W)
          IF (JA.EQ.1) THEN
            DO 1380 I = 1 , N
              T1 = T0
              CALL DRKFGS(FDIF,N,WF(1,I),T1,HS,5,W)
 1380       CONTINUE
          ENDIF
            T0 = WGR(MWGR)
            MWGR = MWGR + 1
            IF (T0.NE.T2) GOTO 1340
        ENDIF
C     IF JA=1 COMPUTATION OF JK, QK AND UK
        IF (JA.EQ.1) THEN
          DO 1400 I = 1 , N
          DO 1400 J = 1 , N
            WF(J,I) = (WF(J,I) - S(J,K)) / EPS
            IF (K.EQ.2) WF0(J,I) = WF(J,I)
 1400     CONTINUE
          CALL DQUDEC(WF,N,N,N,N,W(1,1),W(1,2),W(1,3))
          IF (K.EQ.2) THEN
            ISORT = 0
            DO 1700 ITEL = 1 , N
              CALL DSORTD(W(1,1),N,N,KP,IP,ISORT,IER)
              IF (IER.EQ.0) GOTO 1800
              DO 1500 I = 1 , N
                IK = IP(I)
              DO 1500 J = 1 , N
                WF(J,I) = WF0(J,IK)
                Q(J,I,1) = Q(J,IK,NTI)
 1500         CONTINUE
              DO 1600 I = 1 , N
              DO 1600 J = 1 , N
                WF0(I,J) = WF(I,J)
                 Q(I,J,NTI) = Q(I,J,1)
 1600         CONTINUE
              CALL DQUDEC(WF,N,N,N,N,W(1,1),W(1,2),W(1,3))
 1700       CONTINUE
          ENDIF
 1800     CALL DFQUS(WF,N,N,N,W(1,1),W(1,2),K,KKK,U,NU,NTI,Q,WF0,W(1,3))
        ENDIF
 1900 CONTINUE
      RETURN
 5000 IERROR = 122
      RETURN
 5100 IERROR = 210 + IER
      RETURN
 5200 IERROR = 219
      IP(1) = (TI(NRTI) - TI(1)) / (T1 - TI(1)) * (LWGR+1) + 1
      RETURN
C     END OF DCSAOJ
      END
      SUBROUTINE DJINGX(G,N,NRTI,X,S,NTI,Q,U,NU,ER,KKK,JAC,D,KPART,
     1                  PHIREC,W,WF,WF0,WBMA,WBMB,IP1,IP2,IERROR)
C     ---------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N,NTI),S(N,NTI),Q(N,N,NTI),U(NU,NTI),ER(5),D(N,NTI),
     1          PHIREC(NU,NTI),W(N,3),WF(N,N),WF0(N,N),WBMA(N,N),
     2          WBMB(N,N)
      INTEGER KKK(N),IP1(N),IP2(N)
      EXTERNAL G
C
      IERROR = 0
C     COMPUTATION OF THE INHOMOGENEOUS TERMS OF THE RECURSION.
      DO 1400 I = 2 , NRTI
        DO 1100 J = 1 , N
           W(J,1) = S(J,I) - X(J,I)
 1100   CONTINUE
        DO 1300 J = 1 , N
          SOM = 0.D0
          DO 1200 K = 1 , N
            SOM = SOM + Q(K,J,I) * W(K,1)
 1200     CONTINUE
          D(J,I) = SOM
 1300   CONTINUE
 1400 CONTINUE
      IF (JAC.NE.0) THEN
      IP1(1) = 0
      CALL DGTUR(N,1,NRTI,U,NU,NTI,Q,D,ER,KKK,KPART,IP1,IP2,W(1,1),
     1           W(1,2),W(1,3),WF,WF0,IERROR)
      CALL DKPCH(N,U,NU,NTI,NRTI,KPART,KKK,ER,IERROR)
      CALL DFUNRC(N,KPART,NRTI,U,NU,NTI,KKK,PHIREC,W(1,1),W(1,2),IERROR)
      IF (IERROR.NE.0) RETURN
      ENDIF
      CALL DPSR(N,KPART,NRTI,U,NU,NTI,D,KKK,S,W(1,1),W(1,2),IERROR)
      IF (IERROR.NE.0) RETURN
C     COMPUTATION OF THE BOUNDARY CONDITIONS.
      IF (JAC.EQ.0) THEN
        CALL G(N,X(1,1),X(1,NRTI),W(1,1),WF,WF0)
      ELSE
        CALL G(N,X(1,1),X(1,NRTI),W(1,1),WBMA,WBMB)
      ENDIF
      DO 1500 I = 1 , N
        W(I,1) = -W(I,1)
 1500 CONTINUE
      CALL DSBVP(N,1,WBMA,WBMB,NRTI,Q,NTI,W(1,1),PHIREC,NU,S,S,
     1           ER(3),ER(4),KKK,IP1,WF,WF0,W(1,2),W(1,3),IERROR)
      RETURN
C     END OF DJINGX
      END
      SUBROUTINE DNEWPO(FDIF,N,TOL,EPSMA,TI,NTI,NRTI,ALINC,X,S,KKK,Q,U,
     1                  NU,IPG,WG,LWG,W,WF,WF0,IERROR)
C     -----------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION TI(NTI),ALINC(4),X(N,NTI),S(N,NTI),Q(N,N,NTI),
     1          U(NU,NTI),WG(LWG),W(N,7),WF(N,N),WF0(N,N)
      INTEGER KKK(N),IPG(NTI)
      DIMENSION HS(10)
      EXTERNAL FDIF
C
      IERROR = 0
      IF (IPG(1).EQ.0) RETURN
C     INPUT ERROR 122
      IF (NTI-NRTI.LT.1) GOTO 5000
      LHI = 0
      NF = 0
      JAC = 0
      IF (IPG(1).EQ.-1) JAC = 1
      B = TI(NRTI)
      EPS = DSQRT(EPSMA)
      MWGR = 1
      K = 2
C     COMPUTING TOTAL NUMBER OF CURRENT GRID POINTS AND EXPECTED NEEDED
C     NUMBER OF GRIDPOINTS
      NEG = 0
      IJ = 0
      DO 1300 I = 2 , NRTI
        IM1 = I - 1
        DO 1100 J = IJ+1 , LW
          IF (WG(J).EQ.TI(IM1)) IK = J
          IF (WG(J).EQ.TI(I)) GOTO 1200
 1100   CONTINUE
 1200   IJ = J
        IF (IPG(I).NE.0) NEG = NEG + IJ - IK
 1300 CONTINUE
      NTOTGR = IJ
      NEG = NEG + NTOTGR
 1500 IF (K.GE.NTI) GOTO 5000
      NSHPO = 0
      KM1 = K - 1
      KPL1 = K + 1
      T1 = TI(KM1)
      T2 = TI(K)
      IF (T1.GE.B) GOTO 4900
      IF (IPG(K).EQ.0) THEN
        DO 1600 I = 1 , N
          S(I,K) = X(I,KM1)
          IF (JAC.NE.0) THEN
            DO 1550 J = 1 , N
              WF(J,I) = X(J,KM1) + EPS * Q(J,I,KM1)
 1550       CONTINUE
          ENDIF
 1600   CONTINUE
C     AS LONG AS THE IPG(K) ARE ZERO NO JACOBIAN UPDATE IS NECESSARY.
        T0 = T1
 1700   ST = (WG(MWGR) - T0) / 5.D0
        DO 1800 I = 1 , 5
          HS(I) = ST
 1800   CONTINUE
        CALL DRKFGS(FDIF,N,S(1,K),T1,HS,5,W)
        IF (JAC.NE.0) THEN
          DO 1900 I = 1 , N
            T1 = T0
            CALL DRKFGS(FDIF,N,WF(1,I),T1,HS,5,W)
 1900     CONTINUE
        ENDIF
        T0 = WG(MWGR)
        T1 = T0
        MWGR = MWGR + 1
        IF (T1.NE.T2) GOTO 1700
        IF (JAC.NE.0) THEN
          DO 1950 I = 1 , N
          DO 1950 J = 1 , N
            WF(J,I) = (WF(J,I) - S(J,K)) / EPS
 1950     CONTINUE
          CALL DQUDEC(WF,N,N,N,N,W(1,1),W(1,2),W(1,3))
          CALL DFQUS(WF,N,N,N,W(1,1),W(1,2),K,KKK,U,NU,NTI,Q,WF0,W(1,3))
        ENDIF
        K = K + 1
        GOTO 1500
      ENDIF
C     A NEW SHOOTING POINT IS NEEDED AND FROM NOW ON JACOBIAN UPDATES
C     ARE NECESSARY.
      JAC = 1
      NSHPO = 0
      KWGR = MWGR
      KWG = MWGR
      X1 = 0.D0
      DO 2000 I = 1 , N
        X2 = DABS(X(I,KM1))
        IF (X2.GT.X1) X1 = X2
        S(I,K) = X(I,KM1)
        S(I,KPL1) = X(I,KM1)
        Q(I,1,NTI) = 1.D0
        DO 2000 J = 1 , N
          WF(J,I) = X(J,KM1) + EPS * Q(J,I,KM1)
          WF0(J,I) = WF(J,I)
 2000 CONTINUE
      IWGR = 0
      IWG = 0
 2100 CALL DRKFGG(FDIF,N,T1,WG(MWGR),S(1,K),WF,HS,5,W)
C     CHECK INCREMENT
      X2 = 0.D0
      X3 = 0.D0
      DO 2300 I = 1 , N
        ST = DABS(S(I,K))
        IF (ST.GT.X2) X2 = ST
        X4 = 0.D0
        DO 2200 J = 1 , N
          ST = DABS((WF(J,I)-S(I,K))/Q(I,1,NTI))
          IF (ST.GT.X3) X3 = ST
          IF (ST.GT.X4) X4 = ST
 2200   CONTINUE
        Q(I,2,NTI) = X4
 2300 CONTINUE
      IF ((X2/X1.LT.ALINC(1)).AND.(X3.LT.ALINC(3))) THEN
        IWGR = IWGR + 1
        T1 = WG(MWGR)
        MWGR = MWGR + 1
        IF (T1.EQ.T2) GOTO 3500
        GOTO 2100
      ENDIF
      NSHPO = 1
      IF ((IWGR.EQ.0).AND.(WG(MWGR).EQ.T2)) GOTO 2800
      IF (KWGR.EQ.1) THEN
        T1 = TI(KM1)
      ELSE
        T1 = WG(KWGR-1)
      ENDIF
      DO 2400 I = KWGR , KWGR+IWGR
        CALL DRKFGG(FDIF,N,T1,WG(I),S(1,KPL1),WF0,HS,10,W)
        T1 = WG(I)
 2400 CONTINUE
      X4 = 0.D0
      DO 2600 I = 1 , N
        ST = DABS(S(I,KPL1))
        IF (ST.GT.X4) X4 = ST
        X5 = 0.D0
        DO 2500 J = 1 , N
          ST = DABS((WF0(J,I)-S(J,KPL1))/Q(I,1,NTI))
          IF (ST.GT.X5) X5 = ST
 2500   CONTINUE
        Q(I,2,NTI) = Q(I,2,NTI) / X5
 2600 CONTINUE
      X5 = X2 / X4
      DO 2700 I = 1 , N
        IF (Q(I,2,NTI).GT.X5) X5 = Q(I,2,NTI)
 2700 CONTINUE
      IF ((X5.GE.0.9D0).AND.(X5.LE.1.12D0)) GOTO 3200
C     NEW GRID NEEDED.
 2800 JWGR = IWGR + 1
      IF (NTOTGR+JWGR.GT.LWG) GOTO 5200
      DO 2900 I = NTOTGR , MWGR , -1
        WG(I+JWGR) = WG(I)
 2900 CONTINUE
      J = KWGR - 1
      DO 3000 I = JWGR , 1 , -1
        WG(J+2*I) = WG(J+I)
 3000 CONTINUE
      J = KWGR - 2
      IF (KWGR.EQ.1) THEN
        WG(KWGR) = (WG(KWGR+1) + TI(1)) / 2.D0
        IK = 2
      ELSE
        IK = 1
      ENDIF
      DO 3100 I = IK , JWGR
        I1 = J + 2*I
        WG(I1) = (WG(I1+1) + WG(I1-1)) / 2.D0
 3100 CONTINUE
      MWGR = MWGR + JWGR
 3200 IF (WG(MWGR).NE.T2) THEN
        MWGR = MWGR + 1
        KWGR = MWGR
        IWGR = 0
        X1 = 0.D0
        DO 3400 I = 1 , N
          ST = DABS(S(I,K))
          IF (ST.GT.X1) X1 = ST
          S(I,KPL1) = S(I,K)
          X4 = 0.D0
          DO 3300 J = 1 , N
            ST = DABS(WF(J,I))
            IF (ST.GT.X4) X4 = ST
            WF0(J,I) = WF(J,I)
 3300     CONTINUE
          Q(I,1,NTI) = X4
 3400   CONTINUE
        GOTO 2100
      ELSE
        MWGR = MWGR + 1
      ENDIF
 3500 IF (NSHPO.EQ.0) GOTO 4100
C     DETERMINATION OF THE NEW SHOOTING POINT.
      T1 = TI(KM1)
      T2 = TI(K)
      TP = (T2 + T1) / 2.D0
      E1 = DABS(T2-T1)
      DO 3600 I = KWG , LWG
        ST = DABS(TP-WG(I))
        IF (ST.LT.E1) THEN
          E1 = ST
          JP = I
        ELSE
          GOTO 3700
        ENDIF
 3600 CONTINUE
 3700 TP = WG(JP)
C     INSERT NEW SHOOTING POINT
      NRTI = NRTI + 1
      DO 3800 I = NRTI , KPL1 , -1
        IM1 = I-1
        TI(I) = TI(IM1)
        IPG(I) = IPG(IM1)
        DO 3800 J = 1 , N
          X(J,I) = X(J,IM1)
 3800 CONTINUE
      TI(K) = TP
C     COMPUTATION OF S(.,K) AND WF
      DO 3900 I = 1 , N
        S(I,K) = X(I,KM1)
        DO 3900 J = 1 , N
          WF(J,I) = X(J,KM1) + EPS * Q(J,I,KM1)
 3900 CONTINUE
      DO 4000 I = KWG , JP
        CALL DRKFGG(FDIF,N,T1,WG(I),S(1,K),WF,HS,5,W)
        T1 = WG(I)
 4000 CONTINUE
 4100 DO 4150 I = 1 , N
      DO 4150 J = 1 , N
        WF(J,I) = (WF(J,I) - S(J,K)) / EPS
 4150 CONTINUE
      CALL DQUDEC(WF,N,N,N,N,W(1,1),W(1,2),W(1,3))
      CALL DFQUS(WF,N,N,N,W(1,1),W(1,2),K,KKK,U,NU,NTI,Q,WF0,W(1,3))
      IF (NSHPO.EQ.0) THEN
        K = K + 1
        GOTO 1500
      ENDIF
C     COMPUTATION NEW STARTING VALUE FOR X(.,K).
      X2 = U(1,K)
      DO 4200 I = 2 , N
        IK = KKK(I)
        IF (U(IK,K).LT.X2) X2 = U(IK,K)
 4200 CONTINUE
      ST = DLOG(X2) * (TI(K) - TI(KPL1)) / (TI(K) - TI(KM1))
      ST = DEXP(ST)
      AE = DMIN1(TOL,TOL/ST)
      RE=1.1D-12 + 3.D0 * EPSMA
      T1 = TI(KPL1)
      T2 = TI(K)
      IER = 1
      DO 4300 I = 1 , N
        X(I,K) = X(I,KPL1)
 4300 CONTINUE
      CALL DRKFMS(FDIF,N,X(1,K),T1,T2,RE,AE,EPSMA,WF,NF,HI,LHI,W,IER)
      IF (IER.GT.3) GOTO 5100
      DO 4400 I = 1 , N
        S(I,KPL1) = X(I,K)
        DO 4400 J = 1 , N
          WF(J,I) = X(J,K) + EPS * Q(J,I,K)
 4400 CONTINUE
      T1 = TP
      DO 4500 I = JP+1 , MWGR-1
        CALL DRKFGG(FDIF,N,T1,WG(I),S(1,KPL1),WF,HS,5,W)
        T1 = WG(I)
 4500 CONTINUE
      DO 4600 I = 1 , N
      DO 4600 J = 1 , N
        WF(J,I) = (WF(J,I) - S(J,KPL1)) / EPS
 4600 CONTINUE
      CALL DQUDEC(WF,N,N,N,N,W(1,1),W(1,2),W(1,3))
      CALL DFQUS(WF,N,N,N,W(1,1),W(1,2),KPL1,KKK,U,NU,NTI,Q,WF0,W(1,3))
      K = K + 2
      GOTO 1500
 4900 IF (TI(NRTI).NE.B) GOTO 5000
      RETURN
 5000 IERROR = 122
      RETURN
 5100 IERROR = IER + 210
      RETURN
 5200 IERROR = 219
      IPG(1) = NEG
      RETURN
C     END OF DNEWPO
      END
      SUBROUTINE DRKFMS(F,NEQN,Y,T,TOUT,RE,AE,EPS,PHI,NPHI,
     1                  HI,NHI,W,IFLAG)
C     ----------------------------------------------------------------
C
C     FEHLBERG FOURTH-FIFTH ORDER RUNGE-KUTTA METHOD
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Y(NEQN),PHI(NEQN,NEQN),W(NEQN,7)
      LOGICAL HFAILD,OUTPUT
      EXTERNAL F
C
C     REMIN IS THE MINIMUM ACCEPTABLE VALUE OF RELERR.  ATTEMPTS
C     TO OBTAIN HIGHER ACCURACY WITH THIS SUBROUTINE ARE USUALLY
C     VERY EXPENSIVE AND OFTEN UNSUCCESSFUL.
C
      DATA REMIN/1.D-12/
C
      LHI = 0
      IST = 0
      IUST = 0
      KFLAG = 0
C
C
C     CHECK INPUT PARAMETERS
C
C
      IF (NEQN .LT. 1) GO TO 10
      IF ((RE .LT. 0.0D0)  .OR.  (AE .LT. 0.0D0)) GO TO 10
      MFLAG=IABS(IFLAG)
      IF (MFLAG.NE.1) GOTO 10
C
      U26 = 26.0D0*EPS
      GO TO 50
C
C     INVALID INPUT
   10 IFLAG=8
      RETURN
C
   50 RER=2.0D0*EPS+REMIN
      IF (RE .GE. RER) GO TO 55
C
C     RELATIVE ERROR TOLERANCE TOO SMALL
      RE=RER
      KFLAG=3
C
   55 DT=TOUT-T
C
      A=T
      CALL F(A,Y,W(1,1))
      IF (T .NE. TOUT) GO TO 65
      IFLAG=2
      IF (KFLAG.EQ.3) IFLAG = 3
      RETURN
C
C
   65 IF (IFLAG.EQ.1) THEN
      H=DABS(DT)
      TOLN=0.
      DO 70 K=1,NEQN
        S = RE * DABS(Y(K))
        TOL = S + DMAX1(RE,AE)
        IF (TOL .LE. 0.D0) GO TO 70
        TOLN=TOL
        YPK=DABS(W(K,1))
        IF (YPK*H**5 .GT. TOL) H=(TOL/YPK)**0.2D0
   70 CONTINUE
      IF (TOLN .LE. 0.0D0) H=0.0D0
      H=DMAX1(H,U26*DMAX1(DABS(T),DABS(DT)))
      HI = H
      ELSE
      H = DMIN1(DABS(HI),DABS(DT))
      ENDIF
C
C     SET STEPSIZE FOR INTEGRATION IN THE DIRECTION FROM T TO TOUT
C
      H=DSIGN(H,DT)
C
C     INITIALIZE OUTPUT POINT INDICATOR
C
      OUTPUT= .FALSE.
C
C     TO AVOID PREMATURE UNDERFLOW IN THE ERROR TOLERANCE FUNCTION,
C     SCALE THE ERROR TOLERANCES
C
      SCALE=2.0D0/RE
      AAE=SCALE*AE
C
C
C     STEP BY STEP INTEGRATION
C
  100 HFAILD= .FALSE.
C
C     SET SMALLEST ALLOWABLE STEPSIZE
C
      HMIN=U26*DABS(T)
C
C     ADJUST STEPSIZE IF NECESSARY TO HIT THE OUTPUT POINT.
C     LOOK AHEAD TWO STEPS TO AVOID DRASTIC CHANGES IN THE STEPSIZE AND
C     THUS LESSEN THE IMPACT OF OUTPUT POINTS ON THE CODE.
C
      DT=TOUT-T
      IF (DABS(DT) .GE. 2.0D0*DABS(H)) GO TO 200
      IF (DABS(DT) .GT. DABS(H)) GO TO 150
C
C     THE NEXT SUCCESSFUL STEP WILL COMPLETE THE INTEGRATION TO THE
C     OUTPUT POINT
C
      OUTPUT= .TRUE.
      H=DT
      GO TO 200
C
  150 H=0.5D0*DT
C
C     ADVANCE AN APPROXIMATE SOLUTION OVER ONE STEP OF LENGTH H
C
  200 CALL DRKF1S(F,NEQN,Y,T,H,W)
      IST = IST + 1
C
C     COMPUTE AND TEST ALLOWABLE TOLERANCES VERSUS LOCAL ERROR ESTIMATES
C     AND REMOVE SCALING OF TOLERANCES. NOTE THAT RELATIVE ERROR IS
C     MEASURED WITH RESPECT TO THE AVERAGE OF THE MAGNITUDES OF THE
C     SOLUTION AT THE BEGINNING AND END OF THE STEP.
C
      EEOET=0.0D0
      DO 250 K=1,NEQN
        ET=DABS(Y(K))+DABS(W(K,7))+AAE
        IF (ET .GT. 0.0D0) GO TO 240
C
C       INAPPROPRIATE ERROR TOLERANCE
        IFLAG=5
        RETURN
C
  240   EE=DABS((-2090.0D0*W(K,1)+(21970.0D0*W(K,4)-15048.0D0*W(K,5)))+
     1                        (22528.0D0*W(K,3)-27360.0D0*W(K,6)))
  250   EEOET=DMAX1(EEOET,EE/ET)
C
      ESTTOL=DABS(H)*EEOET*SCALE/752400.0D0
C
      IF (ESTTOL .LE. 1.0D0) GO TO 260
C
C
C     UNSUCCESSFUL STEP
C                       REDUCE THE STEPSIZE , TRY AGAIN
C                       THE DECREASE IS LIMITED TO A FACTOR OF 1/10
C
      IUST = IUST + 1
      HFAILD= .TRUE.
      OUTPUT= .FALSE.
      S=0.1D0
      IF (ESTTOL .LT. 59049.0D0) S=0.9D0/ESTTOL**0.2D0
      H=S*H
      IF (DABS(H) .GT. HMIN) GO TO 200
C
C     REQUESTED ERROR UNATTAINABLE AT SMALLEST ALLOWABLE STEPSIZE
      IFLAG=6
      KFLAG=6
      RETURN
C
C
C     SUCCESSFUL STEP
C                        STORE SOLUTION AT T+H
C                        AND EVALUATE DERIVATIVES THERE
C
  260 LHI = LHI + 1
      DO 270 K = 1 , NEQN
        Y(K) = W(K,7)
  270 CONTINUE
      IF (NPHI.GT.0) THEN
      DO 290 K = 1 , NPHI
        CALL F(T,PHI(1,K),W(1,1))
        CALL DRKF1S(F,NEQN,PHI(1,K),T,H,W)
        DO 280 L = 1 , NEQN
          PHI(L,K) = W(L,7)
  280   CONTINUE
  290 CONTINUE
      ENDIF
      T = T + H
      IF (LHI.EQ.NHI) GOTO 320
      A=T
      CALL F(A,Y,W(1,1))
C
C
C                       CHOOSE NEXT STEPSIZE
C                       THE INCREASE IS LIMITED TO A FACTOR OF 5
C                       IF STEP FAILURE HAS JUST OCCURRED, NEXT
C                          STEPSIZE IS NOT ALLOWED TO INCREASE
C
      S=5.0D0
      IF (ESTTOL .GT. 1.889568D-4) S=0.9D0/ESTTOL**0.2D0
      IF (HFAILD) S=DMIN1(S,1.0D0)
      H=DSIGN(DMAX1(S*DABS(H),HMIN),H)
      HI = H
C
C     END OF CORE INTEGRATOR
C
C     SHOULD WE TAKE ANOTHER STEP
C
      IF (OUTPUT) GO TO 300
      GO TO 100
C
C
C     INTEGRATION SUCCESSFULLY COMPLETED
C
C     INTERVAL MODE
  300 T=TOUT
      IFLAG=2
      IF (KFLAG.EQ.3) IFLAG = 3
      RETURN
C
  320 TOUT = T
      IFLAG = 2
      IF (KFLAG.EQ.3) IFLAG = 3
      RETURN
C
C     END OF DRKFSM
      END
      SUBROUTINE DRKFGG(FDIF,N,T1,T2,X,PHI,HS,NHS,W)
C     ----------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N),PHI(N,N),HS(NHS),W(N,7)
      EXTERNAL FDIF
C
      T0 = T1
      ST = (T2-T1) / NHS
      DO 1000 I = 1 , NHS
        HS(I) = ST
 1000 CONTINUE
      CALL DRKFGS(FDIF,N,X,T1,HS,NHS,W)
      DO 1100 I = 1 , N
        T1 = T0
        CALL DRKFGS(FDIF,N,PHI(1,I),T1,HS,NHS,W)
 1100 CONTINUE
      RETURN
C     END OF DRKFGG
      END
