      SUBROUTINE BNDSCO(F,R,METHOD,X,XS,Y,WORK,IWORK,JS,N,M,MS,
     *      KS,TOL,ITMAX,KP,MMAX,MSMAX,MMS,NMS,NDW,NDIW,NFILE)
C
C
C***********************************************************************
C                                                                      *
C     B O U N D S C O                                                  *
C                                                                      *
C     REALIZATION OF THE MULTIPLE SHOOTING METHOD                      *
C                                                                      *
C     SPECIAL VERSION OF THE PROGRAM  B O U N D S O L                  *
C     FOR THE SOLUTION OF MULTI-POINT BOUNDARY-VALUE PROBLEMS          *
C     WITH SWITCHING CONDITIONS                                        *
C                                                                      *
C     WRITTEN BY   H.J. OBERLE                                         *
C                  UNIVERSITY OF HAMBURG                               *
C                                                                      *
C     REFERENCES  :                                                    *
C                                                                      *
C     BULIRSCH, R.:  DIE MEHRZIELMETHODE ZUR NUMERISCHEN LOESUNG       *
C                    VON NICHTLINEAREN RANDWERTPROBLEMEN UND           *
C                    AUFGABEN DER OPTIMALEN STEUERUNG.                 *
C                    REPORT OF THE CARL CRANZ GESELLSCHAFT,            *
C                    OBERPFAFFENHOFEN, WEST GERMANY, 1971.             *
C     DEUFLHARD, P.: A RELAXATION STRATEGY FOR THE MODIFIED NEWTON     *
C                    METHOD. IN: BULIRSCH, R., OETTLI, W. AND          *
C                    STOER, J. (EDS.) PROCEEDINGS OF CONFERENCE ON     *
C                    OPTIMIZATION AND OPTIMAL CONTROL, OBERWOLFACH.    *
C                    LECTURE NOTES VOL.477, SPRINGER, 1975.            *
C     OBERLE, H.J.:  NUMERISCHE BERECHNUNG OPTIMALER STEUERUNGEN VON   *
C                    HEIZUNG UND KUEHLUNG FUER EIN REALISTISCHES       *
C                    SONNENHAUSMODELL. TECHNISCHE UNIVERSIT[T          *
C                    MUENCHEN, TU-BERICHT NR. M8310, 1983.             *
C     OBERLE, H.J.:  BOUNDSCO - HINWEISE ZUR BENUTZUNG DES MEHRZIEL-   *
C                    VERFAHRENS FUER DIE NUMERISCHE LOESUNG VON        *
C                    RANDWERTPROBLEMEN MIT SCHALTBEDINGUNGEN.          *
C                    HAMBURGER BEITRAEGE ZUR ANGEWANDTEN MATHEMATIK,   *
C                    BERICHTE 6, 1987.                                 *
C                                                                      *
C     UPDATE : SEPTEMBER 1992                                          *
C                                                                      *
C***********************************************************************
C                                                                      *
C --- EXTERNAL SUBROUTINES :                                           *
C                                                                      *
C     F(X,Y,DY,J,L,JS,MMS) ...:                                        *
C                     RIGHT-HAND SIDE OF SYSTEM OF FIRST-ORDER         *
C                     DIFFERENTIAL EQUATIONS                           *
C                     X, Y(N) : ACTUAL INTEGRATION NODE                *
C                     DY(N)   : VALUES OF THE DERIVATIVES AT X         *
C                     J : NUMBER OF SHOOTING INTERVAL (INCLUDING THE   *
C                     SWITCHING POINTS)                                *
C                     L : NUMBER OF SHOOTING INTERVAL (WITHOUT THE     *
C                     SWITCHING POINTS)                                *
C                     JS : INDICATES THE POSITION OF THE               *
C                     ACTUAL INTEGRATION POINT WITH RESPECT TO THE     *
C                     SWITCHING POINTS                                 *
C                     JS(J,K)= 1, IF X.GT.XS(K) , K=1,..,MS            *
C                     JS(J,K)=-1, IF X.LT.XS(K) , K=1,..,MS            *
C                     MMS : DIMENSION PARAMETER                        *
C                                                                      *
C     R(YA,YB,ZZ,W,NYA,NSK,J,L,LS,JS,MMS) ...:                         *
C                     TWO-POINT BOUNDARY CONDITIONS AND                *
C                     SWITCHING CONDITIONS                             *
C                  -  YA, YB: VALUES OF Y(A), Y(B)                     *
C                  -  ZZ : VALUE OF Y(X) AT THE ACTUAL SWITCHING POINT *
C                  -  W(NMS) : BOUNDARY- AND SWITCHING CONDITIONS      *
C                  -  NYA(NMS) : INTEGER ARRAY INDICATING THE STATE    *
C                     VARIABLES FIXED AT THE INITIAL POINT X=A :       *
C                     NYA(I)=K , IF W(K)=YA(I)-PRESCIBED VALUE         *
C                  -  NSK(NMS) : INTEGER ARRAY INDICATING THE          *
C                     SWITCHING CONDITIONS :                           *
C                     NSK(K)=L , IF W(K) DESCRIBES A SWITCHING CONDI-  *
C                     TION, WHICH HAS TO BE SATISFIED AT XS(L).        *
C                  -  LS : INTEGER PARAMETER                           *
C                     LS=0, IF SUBROUTINE R IS CALLED AT X=A, X=B      *
C                     LS=L, IF SUBROUTINE R IS CALLED AT THE SWIT-     *
C                     CHING POINT XS(L), L=1,..MS                      *
C                                                                      *
C     METHOD(N,F,X,Y,XEND,TOL,HMAX,H,J,L,JS,MMS)...:                   *
C                     SUBROUTINE FOR THE SOLUTION OF INITIAL-          *
C                     VALUE PROBLEMS                                   *
C                     METHOD SOLVES THE INITIAL VALUE PROBLEM WITHIN   *
C                     THE INTERVAL (X,XEND) WITH THE (RELATIVE)        *
C                     TOLERANCE TOL.                                   *
C                     H IS THE INITIAL STEPSIZE, AND THE LAST SUCCESS- *
C                     FUL STEPSIZE, RESPECTIVELY.                      *
C                     FAILURE IS INDICATED BY H=0.0                    *
C                                                                      *
C***********************************************************************
C                                                                      *
C --- INPUT PARAMETERS :                                               *
C                                                                      *
C     X(M),Y(N,M) ... INITIAL DATA                                     *
C     XS(MS) ........ INITIAL GUESS OF SWITCHING POINTS                *
C     WORK(NDW) ..... REAL ARRAY USED FOR INTERNAL STORAGE             *
C     IWORK(NDIW) ... INTEGER ARRAY USED FOR INTERNAL STORAGE          *
C     JS(MMS,MSMAX) . INTEGER ARRAY USED FOR THE POSITION OF THE       *
C                     ACTUAL INTEGRATION POINT WITH RESPECT TO         *
C                     SWITCHING POINTS                                 *
C                     JS(J,K)= 1, IF X.GT.XS(K) , K=1,..,MS            *
C                     JS(J,K)=-1, IF X.LT.XS(K) , K=1,..,MS            *
C                     J : SHOOTING INTERVAL (INCLUDING THE             *
C                     SWITCHING POINTS)                                *
C     N ............. NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS     *
C     M ............. NUMBER OF NODES (M.GE.3)                         *
C     MS ............ NUMBER OF SWITCHING POINTS                       *
C     MMAX .......... DIMENSION PARAMETER; M.LE.MMAX                   *
C     MSMAX ......... DIMENSION PARAMETER; MS.LE.MSMAX                 *
C     MMS ........... DIMENSION PARAMETER; M+MS.LE.MMS                 *
C     NMS ........... DIMENSION PARAMETER; N+MS.LE.NMS                 *
C     NDW ........... DIMENSION PARAMETER;                             *
C                       NDW.GE.NMS*(8+MMAX*(10+6*NMS))                 *
C     NDIW .......... DIMENSION PARAMETER; NDIW.GE.5*NMS               *
C     NFILE ......... FILE NUMBER FOR OUTPUT                           *
C     KS ............ 0         NUMBER OF NODES UNCHANGED              *
C                     1         SWITCHING POINTS INSERTED              *
C     TOL ........... REQUIRED RELATIVE PRECISION OF SOLUTION          *
C     ITMAX ......... MAXIMUM PERMITTED NUMBER OF ITERATIONS           *
C                     ITMAX=0   COMPUTATION OF ONE BASIC TRAJECTORY    *
C                               (MULTIPLE SHOOTING)                    *
C                     ITMAX=-1  COMPUTATION OF ONE BASIC TRAJECTORY    *
C                               (SIMPLE SHOOTING)                      *
C     KP ............ PRINT PARAMETER                                  *
C                     KP<-1     NO PRINT                               *
C                     KP=-1     LAST VALUE OF LEVEL FUNCTIONS          *
C                     KP= 0     INITIAL DATA                           *
C                               ITERATIVE VALUES OF LEVEL FUNCTIONS    *
C                               SOLUTION DATA                          *
C                     KP=+1     ADDITIONALLY                           *
C                            ITERATES X(J),Y(I,J),XS(K)                *
C                            ITERATIVE VALUES AT SWITCHING POINTS      *
C                                                                      *
C***********************************************************************
C                                                                      *
C --- OUTPUT PARAMETERS :                                              *
C                                                                      *
C     M ..........    LIKE INPUT M, IF KS=0                            *
C                     NUMBER OF NODES AND SWITCHING POINTS, IF KS=1    *
C     X(M),Y(N,M)     SOLUTION DATA (OR FINAL DATA, RESPECTIVELY)      *
C     XS(MS) .....    SWITCHING POINTS                                 *
C     KP ......... .GT.0     NUMBER OF ITERATIONS PERFORMED            *
C                            TO OBTAIN THE SOLUTION                    *
C                  .LT.0     TEMINATION OF THE ITERATION               *
C                  KP=-1     ITERATION TERMINATES AT STATIONARY POINT  *
C                  KP=-2     ITERATION TERMINATES AFTER ITMAX          *
C                            ITERATIONS                                *
C                  KP=-3     INTEGRATION OF AN INITIAL VALUE PROBLEM   *
C                            TERMINATES                                *
C                  KP=-4     NEWTON METHOD FAILS TO CONVERGE           *
C                  KP=-5     NUMERICAL DIFFERENTIATION IS TERMINATED   *
C                            DUE TO NEIGHBOURING SINGULAR TRAJECTORIES *
C                  KP=-6     SINGULAER ITERATION MATRIX                *
C                  KP=-7     INCONSISTENT ESTIMATION OF THE SWITCHING  *
C                            POINTS                                    *
C                  KP=-8     INCONSISTENT ESTIMATION OF THE INITIAL    *
C                            DATA                                      *
C                                                                      *
C***********************************************************************
C                                                                      *
C --- INTERNAL PARAMETERS :                                            *
C                                                                      *
C     EPMACH ..... RELATIVE MACHINE PRECISION                          *
C     EPMIN ...... MINIMUM PRECISION PARAMETER FOR THE IVP SOLVER      *
C     ETA ........ RELATIVE DEVIATION FOR NUMERICAL DIFFERENTIATION    *
C     FCMIN ...... MINIMUM VALUE FOR THE RELAXATION FACTOR             *
C     COND ....... MAXIMUM LOWER BOUND OF CONDITION NUMBER FOR THE     *
C                  MULTIPLE SHOOTING MATRIX                            *
C     KBROY ...... PARAMETER FOR BROYDEN APPROXIMATIONS                *
C                  0 : ITERATION WITHOUT BROYDEN UPDATING              *
C                  1 : BROYDEN UPDATING                                *
C                                                                      *
C                                                                      *
C***********************************************************************
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      EXTERNAL F,R,METHOD,EXIT,TRAJEC,JACOBI,BROYDN,HOUSE,SUBST
C
      DIMENSION X(*),XS(*),Y(NMS,*),WORK(NDW)
      INTEGER JS(MMS,*), IWORK(NDIW)
C
      DATA EPMACH /1.D-14/, EPMIN /1.D-10/, ETA /1.D-7/,
     *     FCMIN /1.D-2/  , COND /1.D+20/ , HINIT /0.01D0/,
     *     KBROY /1/
C
C --- TEST FOR INVALID INPUT :
C
      IFAIL = 0
C
      IF(MS.GT.MSMAX.OR.M.GT.MMAX.OR.
     *   (M+MS).GT.MMS.OR.(N+MS).GT.NMS.OR.
     *   NDW.LT.NMS*(8+MMAX*(10+6*NMS)).OR.
     *   NDIW.LT.5*NMS) THEN
         IFAIL=1
         GOTO 900
      ENDIF
      IF(M.LT.3) THEN
         IFAIL=4
         GOTO 900
      ENDIF
C
      XX1=X(1)
      XXM=X(M)
      H1=XXM-XX1
C
      DO 100 K=1,M
         T=X(K)
         K1=K+1
         IF(K1.GT.M) K1=K1-2
         H=(X(K1)-X(K)) * 1.D-5
         IF(K.LT.M.AND.H*H1.LE.0.D0) THEN
            IFAIL=2
            GOTO 900
         ENDIF
         IF(MS.GT.0) THEN
         DO 101 I=1,MS
            S=XS(I)
            IF((S-XX1)*(S-XXM).GT.0.D0) THEN
               IFAIL=3
               GOTO 900
            ENDIF
            IF(S.EQ.T) XS(I)=S+H
  101    CONTINUE
         ENDIF
  100 CONTINUE
C
C --- ALLOCATE STORAGE FROM WORK AND IWORK ARRAYS :
C
      IK=NMS*NMS*MMAX
      K1=1
      K2=K1+IK
      K3=K2+IK
      K4=K3+IK
      K5=K4+IK
      K6=K5+IK
      K7=K6+IK
C
      IK=NMS*MMAX
      K8=K7+IK
      K9=K8+IK
      K10=K9+IK
      K11=K10+IK
      K12=K11+IK
      K13=K12+IK
      K14=K13+IK
      K15=K14+IK
      K16=K15+IK
      K17=K16+IK
C
      IK=NMS
      K18=K17+IK
      K19=K18+IK
      K20=K19+IK
      K21=K20+IK
      K22=K21+IK
      K23=K22+IK
      K24=K23+IK
C
      I1=1
      I2=I1+IK
      I3=I2+IK
      I4=I3+IK
      I5=I4+IK
C
      CALL NEWTON (F,R,METHOD,X,XS,Y,TOL,
     *     WORK(K1),WORK(K2),WORK(K3),WORK(K4),WORK(K5),
     *     WORK(K6),WORK(K7),WORK(K8),WORK(K9),WORK(K10),
     *     WORK(K11),WORK(K12),WORK(K13),WORK(K14),WORK(K15),
     *     WORK(K16),WORK(K17),WORK(K18),WORK(K19),WORK(K20),
     *     WORK(K21),WORK(K22),WORK(K23),WORK(K24),
     *     COND,EPMACH,EPMIN,ETA,FCMIN,HINIT,IWORK(I1),
     *     IWORK(I2),IWORK(I3),IWORK(I4),IWORK(I5),JS,IFAIL,ITER,
     *     KBROY,N,M,MS,KS,ITMAX,KP,MMAX,MMS,NMS,NFILE)
C
      GOTO 910
C
C --- INVALID ESTIMATION OF THE INITIAL DATA
C
  900 CONTINUE
      CALL EXIT1(X,XS,Y,KP,IFAIL,N,M,MS,MSMAX,MMAX,MMS,NMS,
     *           NDW,NDIW,NFILE)
  910 CONTINUE
C
C
C***********************************************************************
C     END SUBROUTINE BNDSCO
C
      RETURN
      END




      SUBROUTINE NEWTON (F,R,METHOD,X,XS,Y,TOL,
     *           G,A,GG,AG,HG,SG,DG,DY,DY1,DDY,DY1A,HH,HHA,
     *           YA,YU,YW,U,UH,V,VV,W,WH,Y1,Y2,
     *           COND,EPMACH,EPMIN,ETA,FCMIN,HINIT,NYA,NSK,NSB,
     *           IROW,ICOL,JS,IFAIL,ITER,KBROY,N,M,MS,
     *           KS,ITMAX,KP,MMAX,MMS,NMS,NFILE)
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      EXTERNAL F,R,METHOD,EXIT
C
      DIMENSION X(*),XS(*),Y(NMS,*)
     *,G(NMS,NMS,MMAX),A(NMS,NMS,MMAX),GG(NMS,NMS,MMAX)
     *,AG(NMS,NMS,MMAX),HG(NMS,NMS,MMAX),SG(NMS,NMS,MMAX)
     *,DG(NMS,MMAX),DY(NMS,MMAX),DY1(NMS,MMAX),DDY(NMS,MMAX)
     *,DY1A(NMS,MMAX),HH(NMS,MMAX),HHA(NMS,MMAX),YA(NMS,MMAX)
     *,YU(NMS,MMAX),YW(NMS,MMAX)
     *,U(NMS),UH(NMS),V(NMS),VV(NMS),W(NMS),WH(NMS)
     *,Y1(NMS),Y2(NMS)
      INTEGER JS(MMS,*)
     *,NYA(NMS),NSK(NMS),NSB(NMS),IROW(NMS),ICOL(NMS)
C
C --- INITIAL PREPARATIONS :
C
      IF (IFAIL.GT.0) GOTO 910
C
      FC=FCMIN
C
      IF(TOL.LT.EPMACH*1.D+4) TOL=DABS(EPMACH*1.D+4)
C
      EPINT=TOL * 1.D-02
      IF(EPINT.LT.EPMIN) EPINT=DABS(EPMIN)
C
      CONDR=1.D0/COND
C
      NN=N
      N=N+MS
      M1=M-1
      XX1=X(1)
      XXM=X(M)
C
      IF(MS.GT.0) THEN
      DO 102 I=NN+1,N
      DO 102 K=1,M
         Y(I,K)=XS(I-NN)
  102 CONTINUE
      DO 103 I=1,MMS
      DO 103 K=1,MS
         JS(I,K)=1
  103 CONTINUE
      ENDIF
C
      ITER=0
      IFIN=0
      KT=0
      SUM1A=0.D0
      FCA=FC
      DO 106 I=1,N
         NYA(I)=0
         NSK(I)=0
         NSB(I)=0
         U(I)=1.D0
         Y1(I)=Y(I,1)
         Y2(I)=Y(I,M)
         T=0.D0
         DO 105 K=1,M
            T=T+DABS(Y(I,K))
  105    CONTINUE
         T=T/DBLE(M)
         IF(T.LE.1.D-8) T=1.D0
         DO 106 K=1,M
            YW(I,K)=T
  106 CONTINUE
      LS=-1
      L=0
C
      CALL R(Y1,Y2,U,W,NYA,NSK,J,L,LS,JS,MMS)
C
      HX=(X(2)-X(1))*HINIT
C
      IFAIL=0
      NA=0
      NC=0
      NB=N+1
      DO 113 I=1,N
         IK=NYA(I)
         J=0
  110    CONTINUE
         J=J+1
         IF(I.EQ.NYA(J)) GOTO 111
         IF(J.LT.N) GOTO 110
         NC=NC+1
         IROW(NC)=I
  111    CONTINUE
         IF(IK.NE.0) GOTO 112
         NA=NA+1
         ICOL(NA)=I
         GOTO 113
  112    CONTINUE
         NB=NB-1
         ICOL(NB)=I
         IF(DABS(W(IK)).GT.10.D0*EPMACH) IFAIL=1
  113 CONTINUE
C
C --- OUTPUT OF THE INITIAL TRAJECTORY :
C
      IF(KP.GE.0.OR.IFAIL.EQ.1)
     * CALL EXIT(X,XS,Y,U,YU,FC,SABS,SK,SUM1,SUM2,SUM3,TD,TOL,
     *     X1,ITER,ITMAX,KP,KS,M,MS,NEW,NN,1,IFAIL,NMS,NFILE)
C
C --- SOLUTION OF THE INITIAL VALUE PROBLEMS (BASIC TRAJECTORY) :
C
  200 CONTINUE
C
      CALL TRAJEC
     *   (F,R,METHOD,X,XS,Y,YU,YW,HH,Y1,Y2,W,U,EPINT,H,HC,HX,
     *    SABS,SUMJ,SUMR,X1,JS,IFIN,ISWIT,IS,ITER,ITMAX,JJ,
     *    KONV,KS,M,M1,M11,MS,N,NN,NYA,NSB,NSK,MMS,NMS)
C
      IF(ISWIT.LT.0)
     *CALL EXIT(X,XS,Y,U,YU,FC,SABS,SK,SUM1,SUM2,SUM3,TD,TOL,
     *     X1,ITER,ITMAX,KP,KS,M,MS,NEW,NN,2,IFAIL,NMS,NFILE)
C
      IF(IFIN.EQ.1) GOTO 901
      IF(ISWIT.EQ.-1) THEN
         FC=FC*0.5D0
         GOTO 303
      ENDIF
      IF(ISWIT.EQ.-2) THEN
         IFAIL=3
         GOTO 910
      ENDIF
C
C --- FIRST MONOTONICITY TEST :
C
      SUM1=SUMJ+SUMR
      SABS=SABS+SUMR
      IF(ITER.GT.0) THEN
         KT=1
         IF(SUM1.LE.SUM1A) KTEST=1
         GOTO 810
      ELSE
         SUMJA=SUMJ
         IF(ITMAX.LE.0.OR.NA.EQ.0) GOTO 900
         GOTO 500
      ENDIF
C
C --- UNDERRELAXATION :
C
  300 CONTINUE
C
      IF(KP.GE.0.AND.KT.GT.0)
     *CALL EXIT(X,XS,Y,U,YU,FC,SABS,SK,SUM1,SUM2,SUM3,TD,TOL,
     *     X1,ITER,ITMAX,KP,KS,M,MS,NEW,NN,3,IFAIL,NMS,NFILE)
C
      T=SQRT(SUM3/SUM3A)
      T=(T+FC-1.D0)*2.D0/FC
      T=SQRT(4.D0*T+1.D0)-1.D0
      FC=FC/T
C
  303 CONTINUE
      IF(FC.LT.FCMIN) GOTO 310
      DO 304 K=1,M1
      DO 304 I=1,N
         Y(I,K)=YA(I,K)+FC*DY(I,K)
  304 CONTINUE
      IF(MS.GT.0) THEN
         DO 305 I=1,MS
            XS(I)=Y(NN+I,1)
  305    CONTINUE
         DO 306 I=1,MS
            IF((XS(I)-XX1)*(XS(I)-XXM).GE.0.D0) THEN
               IFAIL=7
               GOTO 910
            ENDIF
  306    CONTINUE
      ENDIF
      GOTO 200
C
C --- FAILTURE OF THE MODIFIED NEWTON METHOD :
C
  310 CONTINUE
      DO 311 I=1,N
         Y(I,M)=YA(I,M)
         IF(I.GT.NN) THEN
            L=I-NN
            XS(L)=YA(I,1)
         ENDIF
         DO 311 K=1,M1
            Y(I,K)=YA(I,K)
            HH(I,K)=HHA(I,K)
  311 CONTINUE
      SABS=SABSA
      SUMR=SUMRA
      KT=-1
      IF(NEW.GT.0.AND.KONV.EQ.0) GOTO 500
      IFAIL=4
      GOTO 910
C
C --- STARTING THE FOLLOWING INTEGRATION STEP :
C
  400 CONTINUE
      IF(KP.GE.0)
     *CALL EXIT(X,XS,Y,U,YU,FC,SABS,SK,SUM1,SUM2,SUM3,TD,TOL,
     *     X1,ITER,ITMAX,KP,KS,M,MS,NEW,NN,4,IFAIL,NMS,NFILE)
C
      KT=0
      IF(FC.LT.FCA.AND.NEW.GT.0) KTEST=0
      FCA=FC
      DO 403 I=1,N
         U(I)=Y(I,M)
         DO 403 K=1,M1
            DY1A(I,K)=DY1(I,K)
            T=(DABS(Y(I,K))+DABS(YA(I,K)))*0.5D0
            IF(T.LT.DABS(DY(I,K))) GOTO 403
            IF(T.LT.ETA) T=1.D0
            YW(I,K)=T
  403 CONTINUE
      IF(KTEST.EQ.3.AND.KBROY.EQ.1) GOTO 600
C
C --- COMPUTATION OF THE JACOBIAN MATRICES G(I,K,JJ) :
C
  500 CONTINUE
      NEW=0
      CALL JACOBI
     *   (F,R,METHOD,G,A,X,XS,Y,YU,YW,HH,U,V,VV,W,WH,Y1,Y2,
     *    EPINT,ETA,H,HX,X1,JS,ICOL,NYA,NSB,NSK,ISWIT,KONV,
     *    M,M1,MS,N,NA,NN,MMS,NMS)
C
      IF(ISWIT.LT.0) THEN
         CALL EXIT(X,XS,Y,U,YU,FC,SABS,SK,SUM1,SUM2,SUM3,TD,TOL,
     *     X1,ITER,ITMAX,KP,KS,M,MS,NEW,NN,5,IFAIL,NMS,NFILE)
         IFAIL=5
         GOTO 910
      ENDIF
      GOTO 700
C
C --- BROYDEN APPROXIMATION OF THE JACOBIAN MATRICES :
C
  600 CONTINUE
      NEW=NEW+1
C
      CALL BROYDN(F,R,G,A,DY,YA,YW,HH,HHA,Y1,U,UH,V,VV,WH,FCA,
     *           JS,M1,N,NYA,NSB,NSK,MMS,NMS)
C
C --- SOLUTION OF THE LINEAR SYSTEM :
C
  700 CONTINUE
      SABSA=SABS
      SUMRA=SUMR
      KTEST=0
      KRED=0
      DO 701 JJ=1,M1
      DO 701 K=1,N
         HHA(K,JJ)=HH(K,JJ)
  701 CONTINUE
  710 CONTINUE
      IF(KT.GT.0) GOTO 810
C
  800 CONTINUE
C
      CALL HOUSE(N,M1,G,A,GG,AG,SG,HG,DG,YW,SK,TD,CONDR,
     *           JF,NMS)
C
      IF(JF.EQ.1) THEN
         IFAIL=6
         GOTO 910
      ENDIF
  810 CONTINUE
C
      CALL SUBST(N,M1,GG,AG,SG,HG,DG,HH,DDY,NMS)
C
      IF(NA.LT.N) THEN
         DO 811 I=NA+1,N
            DDY(ICOL(I),1)=0.D0
  811    CONTINUE
      ENDIF
      SUM3=0.D0
      IF(KT.LE.0.AND.ITER.GT.0) THEN
         IF(KT.EQ.0) SFC1=0.D0
         SFC2=0.D0
      ENDIF
      DO 817 I=1,N
         T=DDY(I,1)
         SUM3=SUM3+T*T
         S=T*YW(I,1)
         IF(KT.GT.0) THEN
            DY1(I,1)=S
            GOTO 816
         ENDIF
         IF(ABS(T).GT.TOL) KONV=0
         IF(ITER.GT.0) THEN
            T=T-DY1A(I,1)/YW(I,1)
            SFC2=SFC2+T*T
            IF(KT.GE.0) THEN
               T=DY(I,1)/YW(I,1)
               SFC1=SFC1+T*T
            ENDIF
         ENDIF
         DY(I,1)=S
         YA(I,1)=Y(I,1)
         YA(I,M)=Y(I,M)
  816    CONTINUE
         U(I)=S
  817 CONTINUE
C
C --- SECOND MONOTONICITY TEST :
C
      IF(KT.GT.0) THEN
         SUM2=SUMJ+SUM3
         IF(SUM2.LE.SUM2A) KTEST=KTEST+1
      ELSE
         SUM2A=SUM3
         IF(ITER.GT.0.AND.KT.EQ.0) SUMJA=0.D0
      ENDIF
C
      IF(M.GT.2) THEN
      DO 830 JI=2,M1
         J1=JI-1
         DO 830 I=1,N
         ST=HH(I,J1)
         T=DDY(I,JI)
         H=YW(I,JI)
         S=T*H
         SUM3=SUM3+T*T
         IF(KT.GT.0) THEN
            DY1(I,JI)=S
            GOTO 830
         ENDIF
         IF(DABS(T).GT.TOL) KONV=0
         IF(ITER.GT.0) THEN
            T=T-DY1A(I,JI)/H
            SFC2=SFC2+T*T
            IF(KT.EQ.0) THEN
               T=ST/H
               SUMJA=SUMJA+T*T
               T=DY(I,JI)/H
               SFC1=SFC1+T*T
            ENDIF
         ENDIF
         DY(I,JI)=S
         YA(I,JI)=Y(I,JI)
  830 CONTINUE
      ENDIF
C
C --- THIRD MONOTONICITY TEST :
C
  840 CONTINUE
      IF(KT.GT.0) THEN
         IF(SUM3.LE.SUM3A) KTEST=KTEST+1
         IF(KTEST) 300,300,400
      ENDIF
C
      SUM1A=SUMR+SUMJA
      SUM2A=SUM2A+SUMJA
      SUM3A=SUM3
      IF(KONV.EQ.1) THEN
         FC=1.D0
         GOTO 850
      ENDIF
      IF(ITER.EQ.0) THEN
         IF(KT.LT.0) FC=FCMIN
         GOTO 850
      ENDIF
      IF(SFC2.LT.SFC1*1.D-6) SFC2=SFC1*1.D-6
      FCS=SQRT(SFC1/SFC2)*FCA
      IF(NEW.GT.0.AND.FCS.LT.FCA*3.D0) GOTO 310
C
C --- COMPUTATION OF THE RELAXATION FACTOR :
C
      FC=1.D0
      IF(FCS.LT.0.7D0) FC=FCS
      IF(FC.GE.FCMIN) GOTO 850
C
      IF(KP.GE.0)
     *CALL EXIT(X,XS,Y,U,YU,FC,SABS,SK,SUM1,SUM2,SUM3,TD,TOL,
     *     X1,ITER,ITMAX,KP,KS,M,MS,NEW,NN,6,IFAIL,NMS,NFILE)
      GOTO 310
C
  850 CONTINUE
      IF(EPINT*SK.GT.1.D0) EPINT=1.D0/SK
      IF(EPINT*TD.GT.1.D0) EPINT=1.D0/TD
      IF(EPINT.LT.EPMIN) EPINT=EPMIN
C
      IF(KP.GE.0)
     *CALL EXIT(X,XS,Y,U,YU,FC,SABS,SK,SUM1A,SUM2A,SUM3A,TD,TOL,
     *     X1,ITER,ITMAX,KP,KS,M,MS,NEW,NN,7,IFAIL,NMS,NFILE)
C
      ITER=ITER+1
      DO 854 JI=1,M1
      DO 854 I=1,N
         Y(I,JI)=YA(I,JI)+FC*DY(I,JI)
  854 CONTINUE
      IF(MS.GT.0) THEN
         DO 855 I=1,MS
            XS(I)=Y(NN+I,1)
  855    CONTINUE
         DO 856 I=1,MS
            IF((XS(I)-XX1)*(XS(I)-XXM).GE.0.D0) THEN
               IFAIL=7
               GOTO 910
            ENDIF
  856    CONTINUE
      ENDIF
      IF(KONV.EQ.1) GOTO 900
      IF(ITER.GT.ITMAX) THEN
         IFAIL=2
         GOTO 910
      ENDIF
      GOTO 200
C
C --- COMPUTATION OF THE LAST TRAJECTORY :
C
  900 CONTINUE
      IFIN=1
      JJ=M11
      H=HC
      GOTO 200
C
C --- EXIT SOLUTION :
C
  901 CONTINUE
      DO 902 I=1,N
         Y(I,M)=U(I)
  902 CONTINUE
      IF(KP.GE.-1)
     *CALL EXIT(X,XS,Y,U,YU,FC,SABS,SK,SUM1A,SUM2A,SUM3A,TD,TOL,
     *     X1,ITER,ITMAX,KP,KS,M,MS,NEW,NN,8,IFAIL,NMS,NFILE)
      KP=ITER
      N=NN
      RETURN
C
C --- EXIT FAILURE :
C
  910 CONTINUE
      IF(KP.GE.0)
     *CALL EXIT(X,XS,Y,U,YU,FC,SABS,SK,SUM1A,SUM2A,SUM3A,TD,TOL,
     *     X1,ITER,ITMAX,KP,KS,M,MS,NEW,NN,9,IFAIL,NMS,NFILE)
      KP=-IFAIL
      N=NN
      RETURN
C
C***********************************************************************
C     END SUBROUTINE NEWTON
C
      END



      SUBROUTINE TRAJEC
     *   (F,R,METHOD,X,XS,Y,YU,YW,HH,Y1,Y2,W,U,EPINT,H,HC,HX,
     *    SABS,SUMJ,SUMR,X1,JS,IFIN,IFAIL,IS,ITER,ITMAX,JJ,
     *    KONV,KS,M,M1,M11,MS,N,NN,NYA,NSB,NSK,MMS,NMS)
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      EXTERNAL F,R,METHOD
C
      DIMENSION X(*),XS(*),Y(NMS,*)
     *,HH(NMS,*),YU(NMS,*),YW(NMS,*)
     *,U(*),VV(18),W(*),Y1(*),Y2(*)
      INTEGER JS(MMS,*)
      INTEGER NYA(*),NSK(*),NSB(*)
C
      IFAIL=0
      IF(IFIN.EQ.1) THEN
         J=JJ
         L=JJ
         GOTO 202
      ENDIF
C
      J=1
      JJ=1
      L=1
      KONV=1
      M11=M1
      IS=0
      SABS=0.D0
      SUMJ=0.D0
      H=HX
      IF(MS.GT.0) THEN
         DO 201 K=1,MS
            JS(1,K)=-1
            IF(X(1).GT.XS(K)) JS(1,K)=1
  201    CONTINUE
      END IF
C
  202 CONTINUE
      J1=J+1
      JJ1=JJ+1
      DO 203 K=1,N
         U(K)=Y(K,JJ)
  203 CONTINUE
      X1=X(JJ)
      X22=X(JJ1)
C
  204 CONTINUE
      X2=X22
      LS=0
      IF(MS.GT.0) THEN
         DO 205 K=1,MS
         X3=XS(K)
         IF((X2-X3)*(X1-X3).GE.0.D0) GOTO 205
         LS=K
         X2=X3
  205    CONTINUE
      END IF
      HMAX=ABS(X2-X1)
      HR=H
C
      CALL METHOD(NN,F,X1,U,X2,EPINT,HMAX,H,J,L,JS,MMS)
C
      CALL F(X2,U,VV,J,L,JS,MMS)
      IF(H.EQ.0.D0) THEN
         IF(ITER.EQ.0) THEN
            IFAIL=-2
         ELSE
            IFAIL=-1
         END IF
         RETURN
      END IF
C
C --- EVALUATION OF THE SWITCHING CONDITIONS :
C
  210 CONTINUE
      IF(MS.EQ.0) GOTO 222
      DO 211 K=1,MS
      JS(J1,K)=JS(J,K)
  211 CONTINUE
      IF(LS.EQ.0) GOTO 222
      IF(IS.NE.0) GOTO 212
      HC=HR
      IS=J
      M11=IS
  212 CONTINUE
C
      CALL R(Y1,Y2,U,W,NYA,NSK,J,L,LS,JS,MMS)
C
      DO 213 K=1,N
         IF(NSK(K).NE.LS) GOTO 213
         HH(K,M1)=W(K)
         NSB(K)=L
  213 CONTINUE
      JS(J1,LS)=-JS(J,LS)
      J=J1
      J1=J+1
      X1=XS(LS)
C
      IF(IFIN*KS.EQ.0) GOTO 204
C
C --- INSERTION OF THE SWITCHING POINTS (LAST TRAJECTORY) :
C
      M=M+1
      MJ=M-J
      JJ=J
      JJ1=J1
      DO 214 K=1,MJ
         JK=M-K
         JI=JK+1
         X(JI)=X(JK)
         DO 214 I=1,N
         Y(I,JI)=Y(I,JK)
  214 CONTINUE
      X(J)=X1
      DO 215 I=1,N
         Y(I,J)=U(I)
  215 CONTINUE
      GOTO 204
C
  222 CONTINUE
      L=L+1
      IF(JJ1.EQ.M) GOTO 240
      IF(IFIN.EQ.1) GOTO 231
C
C --- CONTINUITY CONDITIONS :
C
      DO 230 K=1,N
         T=U(K)
         IF(ITMAX.EQ.-1) Y(K,JJ1)=T
         YU(K,JJ)=T
         T=T-Y(K,JJ1)
         HH(K,JJ)=T
         SABS=SABS+T*T
         T=T/YW(K,JJ1)
         SUMJ=SUMJ+T*T
  230 CONTINUE
  231 CONTINUE
      JJ=JJ1
      J=J1
      GOTO 202
C
C --- BOUNDARY CONDITIONS  :
C
  240 CONTINUE
      DO 241 I=1,N
         Y(I,JJ1)=U(I)
  241 CONTINUE
C
      IF(IFIN.EQ.1) RETURN
C
      IF(M11.EQ.M1) HC=HR
      DO 242 I=1,N
         Y2(I)=U(I)
         Y1(I)=Y(I,1)
  242 CONTINUE
C
      CALL R(Y1,U,U,W,NYA,NSK,J,L,LS,JS,MMS)
C
      SUMR=0.D0
      DO 243 K=1,N
         IF(NSK(K).NE.0) W(K)=HH(K,M1)
         T=W(K)
         HH(K,M1)=T
         SUMR=SUMR+T*T
         Y(K,M)=U(K)
  243 CONTINUE
C
C***********************************************************************
C     END SUBROUTINE TRAJEC
C
      RETURN
      END


      SUBROUTINE JACOBI
     *   (F,R,METHOD,G,A,X,XS,Y,YU,YW,HH,U,V,VV,W,WH,Y1,Y2,
     *    EPINT,ETA,H,HX,X1,JS,ICOL,NYA,NSB,NSK,IFAIL,KONV,
     *    M,M1,MS,N,NA,NN,MMS,NMS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      EXTERNAL F,R,METHOD
C
      DIMENSION X(*),XS(*),Y(NMS,*)
     *,G(NMS,NMS,*),A(NMS,NMS,*)
     *,HH(NMS,*),YU(NMS,*),YW(NMS,*)
     *,U(*),V(*),VV(*),W(*),WH(*),Y1(*),Y2(*)
      INTEGER JS(MMS,*)
      INTEGER NYA(*),NSK(*),NSB(*),ICOL(*)
C
C     COMPUTATION OF THE JACOBIAN MATRIX BY NUMERICAL
C     DIFFERENTIATION
C
      IFAIL=0
      XETA=ETA
C
      DO 501 JJ=1,M
      DO 501 K=1,N
      DO 501 I=1,N
         G(I,K,JJ)=0.D0
         A(I,K,JJ)=0.D0
  501 CONTINUE
C
      J=1
      L=1
      JZ=0
      HR=HX
      HR1=HR
      IF(KONV.NE.1) THEN
         DO 502 I=1,N
            W(I)=HH(I,M1)
            U(I)=Y(I,M)
            Y1(I)=Y(I,1)
  502    CONTINUE
      ENDIF
C
C --- TREATMENT OF SUBINTERVALL NO. J  :
C
  510 J1=J+1
      L1=L+1
      X11=X(L)
      X22=X(L1)
C
      DO 533 IK=1,N
C
      I=ICOL(IK)
C
      IF(J.EQ.1.AND.IK.GT.NA) THEN
         A(NYA(I),I,1)=1.D0
         GOTO 533
      END IF
C
      LL=0
  511 CONTINUE
      H=HR
      X1=X11
      X2=X22
      DO 512 K=1,N
         V(K)=Y(K,L)
  512 CONTINUE
C
C --- VARIATION OF THE SWITCHING POINTS :
C
      IF(I.LE.NN) GOTO 520
      LL=I-NN
      IF(L.NE.M1) G(I,I,L)=1.D0
      SS=XS(LL)
      T=SS*XETA
      IF(ABS(SS).LT.1.D-10) T=YW(I,L)*XETA
      A1=X11-SS
      A2=X22-SS
      IF(A1*A2.GT.0.D0) GOTO 532
      IF(SS.LT.0.D0) T=-T
      IF(T/A1.GT.1.D0) T=0.5D0*A1
      IF(T/A2.GT.1.D0) T=0.5D0*A2
      XS(LL)=T+SS
      V(I)=XS(LL)
      T=1.D0/T
      JZ=0
      GOTO 522
C
C --- VARIATION OF THE COMPONENTS :
C
  520 CONTINUE
      YS=V(I)
      T=YW(I,L)*XETA
      IF(YS.LT.0.D0) T=-T
      ST=YS+T
      V(I)=ST
      T=1.D0/T
      JZ=0
      LS=0
      IF(L.NE.1) GOTO 522
      CALL R(V,U,U,WH,NYA,NSK,J,L,LS,JS,MMS)
      DO 521 K=1,N
         IF(NSK(K).EQ.0)  THEN
         A(K,I,1)=(WH(K)-W(K))*T
         ENDIF
  521 CONTINUE
  522 CONTINUE
      LS=0
      IF(MS.GT.0) THEN
         DO 523 K=1,MS
            X3=XS(K)
            IF((X2-X3)*(X1-X3).LT.0.D0) THEN
               LS=K
               X2=X3
            ENDIF
  523    CONTINUE
      ENDIF
      HMAX=ABS(X2-X1)
C
      CALL METHOD(NN,F,X1,V,X2,EPINT,HMAX,H,J,L,JS,MMS)
C
      IF(H.EQ.0.D0) GOTO 527
      CALL F(X2,V,VV,J,L,JS,MMS)
      HR1=H
  525 CONTINUE
      IF(LS.EQ.0) GOTO 528
C
C --- SWITCHING CONDITIONS :
C
      CALL R(Y1,Y2,V,WH,NYA,NSK,J,L,LS,JS,MMS)
C
      DO 526 K=1,N
         IF(NSK(K).EQ.LS.AND.NSB(K).EQ.L)  THEN
         A(K,I,L)=(WH(K)-W(K))*T
         ENDIF
  526 CONTINUE
      J=J1
      J1=J+1
      JZ=JZ+1
      X1=XS(LS)
      X2=X22
      GOTO 522
C
C  ---INTEGRATION FAILED
C     RECOVERY ATTEMPT BY VARIATION IN THE OPPOSITE DIRECTION :
C
  527 CONTINUE
      IF(XETA.GT.0.D0) THEN
         XETA=-XETA
         J=J-JZ
         J1=J1-JZ
         GOTO 511
      ELSE
         IFAIL=-1
         RETURN
      END IF
C
  528 CONTINUE
      IF(XETA.LT.0.D0) XETA=-XETA
      J=J-JZ
      J1=J1-JZ
      IF(L.NE.M1) GOTO 530
C
      L=M
      CALL R(Y1,V,V,WH,NYA,NSK,J,L,LS,JS,MMS)
      L=M1
      DO 529 K=1,N
         G(K,I,M1)=(WH(K)-W(K))*T
         IF(NSK(K).NE.0) G(K,I,M1)=A(K,I,M1)
  529 CONTINUE
      GOTO 532
C
  530 CONTINUE
      DO 531 K=1,NN
         G(K,I,L)=T*(V(K)-YU(K,L))
  531 CONTINUE
C
  532 CONTINUE
      IF(LL.NE.0) XS(LL)=SS
  533 CONTINUE
      HR=HR1
      J=J1+JZ
      L=L1
      IF(L.NE.M) GOTO 510
C
C***********************************************************************
C     END SUBROUTINE JACOBI
C
      RETURN
      END


      SUBROUTINE BROYDN
     *   (F,R,G,A,DY,YA,YW,HH,HHA,Y1,U,UH,V,VV,WH,FCA,
     *    JS,M1,N,NYA,NSB,NSK,MMS,NMS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      EXTERNAL F,R
C
      DIMENSION G(NMS,NMS,*),A(NMS,NMS,*)
     *,DY(NMS,*),HH(NMS,*),HHA(NMS,*),YA(NMS,*),YW(NMS,*)
     *,U(*),UH(*),V(*),VV(*),WH(*),Y1(*)
      INTEGER JS(MMS,*)
      INTEGER NYA(*),NSK(*),NSB(*)
C
C     BROYDEN APPROXIMATION OF THE JACOBIAN MATRIX
C
      LS=0
      L=0
      H=FCA-1.D0
      DO 601 I=1,N
         V(I)=0.D0
         Y1(I)=YA(I,1)
  601 CONTINUE
      CALL R(Y1,U,U,WH,NYA,NSK,J,L,LS,JS,MMS)
      DO 602 I=1,N
         Y1(I)=HH(I,M1)
         IF(NSK(I).NE.0) WH(I)=Y1(I)
  602 CONTINUE
C
      DO 611 JI=1,M1
         NNB=1
         ST=0.D0
         DO 603 I=1,N
            T=DY(I,JI)/YW(I,JI)
            UH(I)=T
            ST=ST+T*T
            S=HH(I,JI)
            IF(JI.EQ.M1) S=Y1(I)
            U(I)=S+H*HHA(I,JI)
  603    CONTINUE
         ST=ST*FCA
         IF(JI.EQ.1) GOTO 605
         IF(JI.EQ.M1) GOTO 609
         DO 604 K=1,N
            IF(NSB(K).EQ.JI) GOTO 605
  604    CONTINUE
         NNB=0
         GOTO 609
  605    CONTINUE
         DO 607 I=1,N
            S=0.D0
            DO 606 K=1,N
               S=S+A(I,K,JI)*DY(K,JI)
  606       CONTINUE
            V(I)=V(I)+S
            IF(NSB(I).EQ.JI) WH(I)=HHA(I,M1)
  607    CONTINUE
         DO 608 I=1,N
            SA=WH(I)+FCA*V(I)
            VV(I)=Y1(I)-SA
            Y1(I)=SA
  608    CONTINUE
  609    CONTINUE
         DO 610 K=1,N
            T=UH(K)/(ST*YW(K,JI))
            DO 610 I=1,N
            S=U(I)
            IF((NNB.GT.0).AND.(JI.LT.M1))
     *      A(I,K,JI)=A(I,K,JI)+T*VV(I)
            G(I,K,JI)=G(I,K,JI)+S*T
  610    CONTINUE
  611 CONTINUE
      DO 612 I=1,N
         Y1(I)=YA(I,1)
  612 CONTINUE
C
C***********************************************************************
C     END SUBROUTINE BROYD
C
      RETURN
      END



      SUBROUTINE HOUSE(N,M1,G,A,GG,AG,SG,HG,DG,YW,SK,TD,CONDR,
     *           JF,NMS)
C
C     HOUSEHOLDER TRIANGULATION OF THE MULTIPE-SHOOTING MATRIX
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION G(NMS,NMS,*),A(NMS,NMS,*)
      DIMENSION GG(NMS,NMS,*),AG(NMS,NMS,*)
      DIMENSION SG(NMS,NMS,*),HG(NMS,NMS,*)
      DIMENSION DG(NMS,*),YW(NMS,*)
C
      LOGICAL BJ
C
      N1=N-1
      M2=M1-1
C
      DO 10 JJ=1,M1
      DO 10 K=1,N
         S=YW(K,JJ)
         DO 20 I=1,N
            GG(I,K,JJ)=G(I,K,JJ)*S
            AG(I,K,JJ)=A(I,K,JJ)*S
            HG(I,K,JJ)=0.D0
   20    CONTINUE
         JJ1=JJ-1
         IF(JJ.GT.1) HG(K,K,JJ1)=-S
   10 CONTINUE
C
      IF(M2.EQ.0) GOTO 510
C
      DO 500 J=1,M2
C
      J1=J+1
      BJ=J.EQ.M2
C
      DO 100 I=1,N
      DO 100 K=1,N
         SG(I,K,J)=0.D0
  100 CONTINUE
C
      DO 500 I=1,N
C
      I1=I-1
      I11=I+1
C
      SZ=0.D0
      DO 210 K=I,N
         SZ=SZ+GG(I,K,J)*GG(I,K,J)
  210 CONTINUE
      DO 220 K=1,I
         SZ=SZ+HG(I,K,J)*HG(I,K,J)
  220 CONTINUE
      V=SQRT(SZ)
      VG=GG(I,I,J)
      IF(VG.LT.0.D0) V=-V
      DG(I,J)=-V
      C=SZ+VG*V
      IF(C.LT.CONDR) GOTO 900
      GG(I,I,J)=VG+V
      IF(I.EQ.N) GOTO 280
      DO 270 K=I11,N
         Y=0.D0
         DO 230 L=I,N
            Y=Y+GG(I,L,J)*GG(K,L,J)
  230    CONTINUE
         IF(I1.GT.0) THEN
         DO 240 L=1,I1
            Y=Y+HG(I,L,J)*HG(K,L,J)
  240    CONTINUE
         ENDIF
         Y=Y/C
         DO 260 L=I,N
            GG(K,L,J)=GG(K,L,J)-GG(I,L,J)*Y
  260    CONTINUE
         DO 270 L=1,I
            HG(K,L,J)=HG(K,L,J)-HG(I,L,J)*Y
  270 CONTINUE
C
  280 CONTINUE
      IF(BJ) GOTO 400
C
      DO 330 K=1,N
         Y=0.D0
         DO 300 L=I,N
            Y=Y+GG(I,L,J)*SG(K,L,J)
  300    CONTINUE
         DO 310 L=1,I
            Y=Y+HG(I,L,J)*GG(K,L,J1)
  310    CONTINUE
         Y=Y/C
         DO 320 L=I,N
            SG(K,L,J)=SG(K,L,J)-GG(I,L,J)*Y
  320    CONTINUE
         DO 330 L=1,I
            GG(K,L,J1)=GG(K,L,J1)-HG(I,L,J)*Y
  330 CONTINUE
C
  400 CONTINUE
C
      DO 500 K=1,N
         Y=0.D0
         DO 410 L=I,N
            Y=Y+GG(I,L,J)*AG(K,L,J)
  410    CONTINUE
         IF(BJ) GOTO 430
C
         DO 420 L=1,I
            Y=Y+HG(I,L,J)*AG(K,L,J1)
  420    CONTINUE
         GOTO 450
C
  430    CONTINUE
         DO 440 L=1,I
            Y=Y+HG(I,L,J)*GG(K,L,J1)
  440    CONTINUE
C
  450    CONTINUE
         Y=Y/C
         DO 460 L=I,N
            AG(K,L,J)=AG(K,L,J)-GG(I,L,J)*Y
  460    CONTINUE
         IF(BJ) GOTO 480
C
         DO 470 L=1,I
            AG(K,L,J1)=AG(K,L,J1)-HG(I,L,J)*Y
  470    CONTINUE
         GOTO 500
C
  480    CONTINUE
         DO 490 L=1,I
            GG(K,L,J1)=GG(K,L,J1)-HG(I,L,J)*Y
  490    CONTINUE
C
  500 CONTINUE
C
  510 CONTINUE
      IF(N.GT.1) THEN
C
      DO 540 I=1,N1
         I11=I+1
         SZ=0.D0
         DO 520 L=I,N
            SZ=SZ+GG(I,L,M1)*GG(I,L,M1)
  520    CONTINUE
         V=SQRT(SZ)
         VG=GG(I,I,M1)
         IF(VG.LT.0.D0) V=-V
         DG(I,M1)=-V
         C=SZ+V*VG
         IF(C.LT.CONDR) GOTO 900
         GG(I,I,M1)=VG+V
         DO 540 K=I11,N
         Y=0.D0
         DO 530 L=I,N
            Y=Y+GG(I,L,M1)*GG(K,L,M1)
  530    CONTINUE
         Y=Y/C
         DO 540 L=I,N
            GG(K,L,M1)=GG(K,L,M1)-GG(I,L,M1)*Y
  540 CONTINUE
      ENDIF
      DG(N,M1)=GG(N,N,M1)
C
C     LOWER ESTIMATION OF NORM AND CONDITION NUMBER
C
      TD=DABS(DG(1,1))
      SK=TD
      DO 600 J=1,M1
      DO 600 I=1,N
         V=ABS(DG(I,J))
         IF(V.LT.SK) SK=V
         IF(V.GT.TD) TD=V
  600 CONTINUE
      IF(SK.LT.CONDR) GOTO 900
      SK=TD/SK
      JF=0
      RETURN
C
  900 CONTINUE
      JF=1
C
C***********************************************************************
C     END SUBROUTINE HOUSE
C
      RETURN
      END



      SUBROUTINE SUBST(N,M1,GG,AG,SG,HG,DG,B,X,NMS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION GG(NMS,NMS,*),AG(NMS,NMS,*)
      DIMENSION SG(NMS,NMS,*),HG(NMS,NMS,*)
      DIMENSION DG(NMS,*),X(NMS,*),B(NMS,*)
C
      LOGICAL BN
C
      BN=N.EQ.1
      M2=M1-1
C
      X(1,1)=-B(1,1)/DG(1,1)
      IF(BN) GOTO 200
      DO 110 I=2,N
         I1=I-1
         V=-B(I,1)
         DO 100 L=1,I1
            V=V-GG(I,L,1)*X(L,1)
  100    CONTINUE
         X(I,1)=V/DG(I,1)
  110 CONTINUE
C
  200 CONTINUE
      IF(M2.LE.0) GOTO 300
      DO 240 J=2,M2
         J1=J-1
         V=-B(1,J)
         DO 210 L=1,N
            V=V-SG(1,L,J1)*X(L,J1)
  210    CONTINUE
         X(1,J)=V/DG(1,J)
         IF(.NOT.BN) THEN
         DO 230 I=2,N
            V=-B(I,J)
            DO 220 L=1,N
               IF(L.LT.I) V=V-GG(I,L,J)*X(L,J)
               V=V-SG(I,L,J1)*X(L,J1)
  220       CONTINUE
            X(I,J)=V/DG(I,J)
  230    CONTINUE
         ENDIF
  240 CONTINUE
C
      DO 260 I=1,N
         V=-B(I,M1)
         DO 250 J=1,M2
         DO 250 L=1,N
            V=V-AG(I,L,J)*X(L,J)
  250    CONTINUE
         X(I,M1)=V
  260 CONTINUE
C
  300 CONTINUE
      X(1,M1)=X(1,M1)/DG(1,M1)
      IF(BN) GOTO 400
      DO 320 I=2,N
         I1=I-1
         V=X(I,M1)
         DO 310 L=1,I1
            V=V-GG(I,L,M1)*X(L,M1)
  310    CONTINUE
         X(I,M1)=V/DG(I,M1)
  320 CONTINUE
C
      N1=N-1
      DO 340 I=N1,1,-1
         C=-DG(I,M1)*GG(I,I,M1)
         Y=0.D0
         DO 330 K=I,N
            Y=Y+GG(I,K,M1)*X(K,M1)
  330   CONTINUE
         Y=Y/C
         DO 340 K=I,N
            X(K,M1)=X(K,M1)-Y*GG(I,K,M1)
  340 CONTINUE
C
  400 CONTINUE
      IF(M2.EQ.0) RETURN
      DO 440 J=M2,1,-1
         J1=J+1
         DO 440 I=N,1,-1
         I1=I-1
         C=-DG(I,J)*GG(I,I,J)
         Y=0.D0
         DO 410 K=1,N
            IF(K.LE.I) Y=Y+HG(I,K,J)*X(K,J1)
            IF(K.GE.I) Y=Y+GG(I,K,J)*X(K,J)
  410    CONTINUE
         Y=Y/C
         DO 420 K=I,N
            X(K,J)=X(K,J)-Y*GG(I,K,J)
  420    CONTINUE
         DO 430 K=1,I
            X(K,J1)=X(K,J1)-Y*HG(I,K,J)
  430    CONTINUE
  440 CONTINUE
C
C***********************************************************************
C     END SUBROUTINE SUBST
C
      RETURN
      END



      SUBROUTINE EXIT
     *  (X,XS,Y,U,YU,FC,SABS,SK,SUM1,SUM2,SUM3,TD,TOL,
     *   X1,ITER,ITMAX,KP,KS,M,MS,NEW,NN,ITEIL,IFAIL,NMS,NFILE)
C
C  ********************************************************************
C                                                                     *
C     SUBROUTINE FOR DATA OUTPUT                                      *
C     EXIT IS CALLED BY THE ROUTINE NEWTON                            *
C                                                                     *
C  ********************************************************************
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION X(*),XS(*),Y(NMS,*),YU(NMS,*),U(*)
C
      GOTO (100,200,300,300,210,400,500,600,700) ITEIL
C
  100 CONTINUE
      WRITE(NFILE,8888)
      WRITE(NFILE,1100)
      DO 110 K=1,M
         WRITE(NFILE,1200) X(K), (Y(I,K),I=1,NN)
  110 CONTINUE
      IF(MS.GT.0) THEN
         WRITE(NFILE,1300)
         WRITE(NFILE,1400) (XS(J),J=1,MS)
      ENDIF
      WRITE(NFILE,1500) NN,M,MS,TOL,ITMAX
      IF(IFAIL.EQ.1) WRITE(NFILE,7100)
      IF(KP.GT.0.OR.ITMAX.LE.0) GOTO 900
      WRITE(NFILE,8888)
      WRITE(NFILE,9999)
      WRITE(NFILE,2100)
      GOTO 900
C
  200 CONTINUE
      WRITE(NFILE,2000) X1,(U(I),I=1,NN)
      GOTO 900
C
  210 CONTINUE
      WRITE(NFILE,2300) X1,(U(I),I=1,NN)
      GOTO 900
C
  300 CONTINUE
C      IF(KP.GT.0) THEN
C         WRITE(NFILE,9999)
C         WRITE(NFILE,2100)
C      ENDIF
      WRITE(NFILE,2200) SABS,SUM1,SUM2,SUM3,FC
      IF(KP.GT.0) WRITE(NFILE,9999)
C
      IF(KP.GT.0.AND.ITEIL.EQ.4) THEN
      DO 310 K=1,M
         WRITE(NFILE,1200) X(K),(Y(I,K),I=1,NN)
  310 CONTINUE
      IF(MS.GT.0) WRITE(NFILE,1400) (XS(K),K=1,MS)
      ENDIF
      GOTO 900
C
  400 CONTINUE
      WRITE(NFILE,4000) ITER,SUM3,FC,SK
      GOTO 900
C
  500 CONTINUE
      IF(KP.GT.0) THEN
         WRITE(NFILE,9999)
         WRITE(NFILE,2100)
      ENDIF
      WRITE(NFILE,5000) ITER,SABS,SUM1,SUM2,SUM3,NEW,SK,TD
      IF(KP.GT.0) WRITE(NFILE,9999)
      GOTO 900
C
  600 CONTINUE
      IF(KP.EQ.-1) THEN
         WRITE(NFILE,8888)
         WRITE(NFILE,9999)
         WRITE(NFILE,2100)
         WRITE(NFILE,5000) ITER,SABS,SUM1,SUM2,SUM3,NEW,SK,TD
         WRITE(NFILE,9999)
         IF(MS.GT.0) THEN
            WRITE(NFILE,1300)
            WRITE(NFILE,1400) (XS(J),J=1,MS)
         ENDIF
         GOTO 900
      ENDIF
      IF(ITMAX.EQ.0) GOTO 603
      IF(KP.EQ.0) WRITE(NFILE,9999)
      WRITE(NFILE,8888)
  603 CONTINUE
      WRITE(NFILE,6100) ITER
      DO 604 J=1,M
         WRITE(NFILE,1200) X(J),(Y(I,J),I=1,NN)
  604 CONTINUE
      IF(MS.GT.0) THEN
         WRITE(NFILE,1300)
         WRITE(NFILE,1400) (XS(J),J=1,MS)
      ENDIF
      WRITE(NFILE,8888)
      WRITE(NFILE,9999)
      IF(ITMAX.NE.0) GOTO 900
C
      WRITE(NFILE,6200) SUM1,SABS
      WRITE(NFILE,9999)
      M2=M-2
      IF(M2.LE.0) GOTO 900
      WRITE(NFILE,6300)
      K=1
      DO 608 JJ=1,M2
         JJ1=JJ+1
         IF(MS.EQ.0.OR.KS.EQ.0) GOTO 607
         DO 606 I=1,MS
            IF(X(JJ1).EQ.XS(I)) GOTO 607
  606    CONTINUE
  607    WRITE(NFILE,1200) X(JJ1),(YU(I,K),I=1,NN)
         K=K+1
  608 CONTINUE
      GOTO 900
C
  700 CONTINUE
      GOTO (720,720,730,740,750,760,770,780) IFAIL
  720 CONTINUE
      WRITE(NFILE,7200) ITMAX
      GOTO 800
  730 CONTINUE
      WRITE(NFILE,7300)
      GOTO 800
  740 CONTINUE
      WRITE(NFILE,7400)
      GOTO 800
  750 CONTINUE
      WRITE(NFILE,7500)
      GOTO 800
  760 CONTINUE
      WRITE(NFILE,7600)
      GOTO 800
  770 CONTINUE
      WRITE(NFILE,7700)
      GOTO 800
  780 CONTINUE
      WRITE(NFILE,7800)
  800 CONTINUE
      WRITE(NFILE,9999)
      DO 810 K=1,M
      WRITE(NFILE,1200) X(K),(Y(I,K),I=1,NN)
  810 CONTINUE
      IF(MS.GT.0) THEN
         WRITE(NFILE,1300)
         WRITE(NFILE,1400) (XS(K),K=1,MS)
      ENDIF
      WRITE(NFILE,8888)
      WRITE(NFILE,9999)
C
  900 CONTINUE
C
C     FORMAT SPECIFICATIONS
C
 1100 FORMAT(' ',78('*')///4X,'Initial Data:'//)
 1200 FORMAT(3X,D15.8,4(D15.8),3(/18X,4(D15.8)))
 1300 FORMAT(/4X,'Switching Points:'/)
 1400 FORMAT(3(/3X,5D15.8))
 1500 FORMAT(//4X,'N=',I2,6X,'M=',I2,6X,'MS=',I2/4X,
     *'Prescribed Relative Precision',D10.2/4X,
     *'Maximum Permitted Number of Iterations',I3///)
 2000 FORMAT(' ******Integration Routine fails at:'/
     *4X,D20.10,5(/24X,4D15.8))
 2100 FORMAT(2X,'It',3X,'Abs.Err.',2X,'Level1',4X,'Level2',4X,'Level3',
     *4X,'Relax. New',2X,'Cond(M)',3X,'Norm(M)')
 2200 FORMAT(6X,D8.2,2X,D8.2,2X,D8.2,2X,D8.2,F7.3)
 2300 FORMAT(' ***** Integration Routine fails in  Numerical Differentiat
     *ion:'/1X,D20.10,4D15.8,5(/21X,4D15.8))
 4000 FORMAT(2X,I2,3X,' Not accepted pre-estimate ',4X,D8.2,F7.3,6X,
     *D8.2)
 5000 FORMAT(2X,I2,2X,D8.2,2X,D8.2,2X,D8.2,2X,D8.2,10X,I2,2X
     *,D8.2,2X,D8.2)
 6100 FORMAT(4X,'Solution obtained after',I3,' iteration step(s)'//,4X,
     *'Solution Data:'/)
 6200 FORMAT(/4X,'Rel. Error',D11.3/4X,'Abs. Error',D11.3)
 6300 FORMAT(4X,'Solution of the Initial Value Problems'/)
 7100 FORMAT(' ******WARNING******:'/
     *       ' Initial data and boundary conditions are inconsistent')
 7200 FORMAT(4X,'Iteration terminates after ITMAX=',I3,' Iterations')
 7300 FORMAT(4X,'BOUNDSCO terminates at a singular point'/4X,'new step
     *size H=0.0')
 7400 FORMAT(4X,'Modified Newton Method fails to converge'/)
 7500 FORMAT(4X,'Numerical Differentiation terminated')
 7600 FORMAT(4X,'Singular Iteration Matrix')
 7700 FORMAT(4X,'Inconsistent estimation of the switching points')
 7800 FORMAT(4X,'Inconsistent estimation of the initial data'/)
 8888 FORMAT(//)
 9999 FORMAT(' ',78('*'))
C
C***********************************************************************
C     END OF SUBROUTINE EXIT
C
      RETURN
      END



      SUBROUTINE EXIT1
     *           (X,XS,Y,KP,IFAIL,N,M,MS,MSMAX,MMAX,MMS,NMS,
     *            NDW,NDIW,NFILE)
C  ********************************************************************
C                                                                      *
C     SUBROUTINE FOR DATA OUTPUT IN CASE OF INVALID INITIAL DATA       *
C     EXIT1 IS CALLED BY THE ROUTINE BNDSCO                            *
C                                                                      *
C  ********************************************************************
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION X(*),XS(*),Y(NMS,*)
C
      GOTO (100,200,300,400) IFAIL
C
  100 CONTINUE
      WRITE(NFILE,8888)
      WRITE(NFILE,1000)
      WRITE(NFILE,1100) N,M,MS,MMAX,MSMAX,MMS,NMS,NDW,NDIW
      GOTO 900
C
  200 CONTINUE
      WRITE(NFILE,8888)
      WRITE(NFILE,2000)
      GOTO 900
C
  300 CONTINUE
      WRITE(NFILE,8888)
      WRITE(NFILE,3000)
C
  400 CONTINUE
      WRITE(NFILE,8888)
      WRITE(NFILE,4000)
C
  900 CONTINUE
      WRITE(NFILE,8888)
      WRITE(NFILE,9100)
      DO 910 K=1,M
         WRITE(NFILE,9200) X(K), (Y(I,K),I=1,N)
  910 CONTINUE
      IF(MS.GT.0) THEN
         WRITE(NFILE,9300)
         WRITE(NFILE,9400) (XS(J),J=1,MS)
      ENDIF
      WRITE(NFILE,8888)
      WRITE(NFILE,9999)
C
      KP=-8
      IF(IFAIL.EQ.3) KP=-7
C
C     FORMAT SPECIFICATIONS
C
 1000 FORMAT(//2X,'Inconsistent estimation of the initial data'/
     *2X,'enlarge dimension of work arrays')
 1100 FORMAT(//2X,'N=',I2,8X,'M=',I2,8X,'MS=',I2/2X,
     *'MMAX=',I2,5X,'MSMAX=',I2,4X,'MMS=',I2/2X,
     *'NMS=',I2,6X,'NDW=',I6,2X,'NDIW=',I6///)
 2000 FORMAT(//2X,'Inconsistent estimation of the initial data'/
     *2X,'The multiple shooting nodes should be ordered in'/
     *2X,'ascending or descending direction')
 3000 FORMAT(2X,'Inconsistent estimation of the switching points')
 4000 FORMAT(2X,'M should be ge 3') 
 9100 FORMAT(' ',78('*')///2X,'Initial Data'//)
 9200 FORMAT(D16.8,4(D16.8),3(/16X,4(D16.8)))
 9300 FORMAT(/2X,'Switching Points'/)
 9400 FORMAT(3(/5D16.8))
 8888 FORMAT(////)
 9999 FORMAT('',78('*'))
C
C  *********************************************************************
C     END OF SUBROUTINE EXIT1
C
      RETURN
      END



      SUBROUTINE DIFSYB(N,F,X,Y,XEND,EPS,HMAX,H,
     *                  JIN,LIN,JSIN,MMS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     EXTRAPOLATION METHOD DUE TO BULIRSCH, GRAGG, STOER
C
C     VERSION: FEBRUARY 1986
C
      DIMENSION Y(N),YA(50),YL(50),YM(50),DY(50),DZ(50),DT(50,10),
     * D(10),S(50),EP(5)
      INTEGER JSIN(MMS,*)
C
      LOGICAL   KONV,BO,KL,GR
      DATA EP /.4D-1,.16D-2,.64D-4,.256D-5,.1024D-6/
      DATA EPMACH /1.D-14/
C
      JMAX=300
      JDIFS=0
      XEPS=(ABS(X)+ABS(XEND))*EPMACH
      H1=XEND-X
      ETA=ABS(EPS)
  401 IF(ABS(H1).LE.XEPS) GO TO 403
      Q=H1/H
      IF(Q.LE.EPMACH) GO TO 403
      IF(Q.GE.1.01D0) GO TO 402
      HR=H
      H=H1
  402 JTI=0
      FY=10.D0
      DO 100 I=1,N
  100 YA(I)=Y(I)
      CALL F(X,Y,DZ,JIN,LIN,JSIN,MMS)
   10 XN=X+H
      BO=.FALSE.
      DO 110 I=1,N
  110 S(I)=0.D0
      M=1
      JR=2
      JS=3
      DO 260 J=1,10
      IF(.NOT.BO) GOTO 200
      D(2)=1.7777777777778D0
      D(4)=7.1111111111111D0
      D(6)=28.444444444444D0
      D(8)=113.77777777778D0
      GOTO 201
  200 D(2)=2.25D0
      D(4)=9.D0
      D(6)=36.D0
      D(8)=144.D0
  201 IF(J.LE.9) GOTO 202
      L=9
      D(9)=256.D0
      GOTO 203
  202 L=J
      D(L)=M*M
  203 KONV=L.GT.3
      M=M+M
      G=H/M
      B=G+G
      DO 210 I=1,N
      YL(I)=YA(I)
  210 YM(I)=YA(I)+G*DZ(I)
      M=M-1
      DO 220 K=1,M
      CALL  F(X+K*G,YM,DY,JIN,LIN,JSIN,MMS)
      DO 220 I=1,N
      U=YL(I)+B*DY(I)
      YL(I)=YM(I)
      YM(I)=U
      U=ABS(U)
      IF(U.GT.S(I)) S(I)=U
  220 CONTINUE
      CALL  F(XN,YM,DY,JIN,LIN,JSIN,MMS)
      KL=L.LT.2
C***  GR=L.GT.5
      GR=L.GT.6
      FS=0.D0
      DO 233 I=1,N
      V=DT(I,1)
      C=(YM(I)+YL(I)+G*DY(I))*0.5D0
      DT(I,1)=C
      TA=C
      IF(KL) GOTO 233
      DO 231 K=2,L
      B1=D(K)*V
      B=B1-C
      W=C-V
      U=V
      IF(B.EQ.0.D0)  GOTO 230
      B=W/B
      U=C*B
      C=B1*B
  230 V=DT(I,K)
      DT(I,K)=U
  231 TA=U+TA
      IF(.NOT.KONV) GOTO 232
      IF (ABS(Y(I)-TA).GT.S(I)*ETA) KONV=.FALSE.
  232 IF(GR.OR.S(I).EQ.0.D0) GOTO 233
      FV=ABS(W)/S(I)
      IF(FS.LT.FV) FS=FV
  233 Y(I)=TA
      IF(FS.EQ.0.D00) GOTO 250
      FA=FY
      K=L-1
      FY=(EP(K)/FS)**(1.D0/(L+K+2.D0*JTI))
      IF(L.EQ.2) GOTO 240
      IF(FY.LT.0.7D0*FA) GOTO 250
  240 IF(FY.GT.0.7D0) GOTO 250
      JTI=JTI+1
      IF(JTI.GT.5) GO TO 250
      H=H*FY
      GOTO 10
  250 IF(KONV) GOTO 20
      D(3)=4.D0
      D(5)=16.D0
      D(7)=64.D0
      BO=.NOT.BO
      M=JR
      JR=JS
  260 JS=M+M
      IF(JTI.GT.5) GOTO 30
      H=H*0.5D0
      GOTO 10
   20 X=XN
      H=H*FY
      H1=XEND-X
      JDIFS=JDIFS+1
      IF(JDIFS.GT.JMAX) GOTO 30
      GOTO 401
  403 H=HR
      RETURN
   30 H=0.D0
      DO 300 I=1,N
  300 Y(I)=YA(I)
      RETURN
      END



      SUBROUTINE RKF7B (N,F,X,Y,XEND,TOL,HMAX,HI,
     *                  JIN,LIN,JSIN,MMS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C     RUNGE-KUTTA SUBROUTINE BASED ON FEHLBERG'S SEVENTH-ORDER FORMULAS
C
      DIMENSION Y(*),
     *          W1(20),W2(20),W3(20),W4(20),W5(20),W6(20),
     *          W7(20),W8(20),W9(20),W10(20),W11(20),W12(20)
      INTEGER JSIN(MMS,*)
      DATA UROUND/1.D-14/
      HMIN=HMAX*TOL*0.01D0
      IND=1
C
C     PREPARE TO MAKE CALCULATIONS
C
100   GOTO (110,120,130), IND
110   CALL F(X,Y,W1,JIN,LIN,JSIN,MMS)
      HMAG=HI
      GOTO 131
120   CALL F(X,Y,W1,JIN,LIN,JSIN,MMS)
      HMAG=1.5D0*ABS(H)
      IF(EST.GT.(0.0279936D0*TOL))
     *   HMAG=0.9D0*(TOL/EST)**(1.D0/7.D0)*ABS(H)
      HMAG=DMAX1(HMAG,HMIN)
      HMAG=DMIN1(HMAG,HMAX)
      GOTO 131
130   HMAG=DMAX1(0.5D0,0.9D0*(TOL/EST)**(1.D0/7.D0))*ABS(H)
      HMAG=DMAX1(HMAG,HMIN)
131   CONTINUE
      IF(HMAG.LT.ABS(XEND-X)) GOTO 140
      GOTO 141
  140 CONTINUE
      H=SIGN(DMIN1(HMAG,0.5D0*ABS(XEND-X)),XEND-X)
      XTRIAL=X+H
      GOTO 142
141   H=XEND-X
      XTRIAL=XEND
142   CONTINUE
C
C     CALCULATIONS
C
      DO 200 I=1,N
200   W12(I)=Y(I)+H*W1(I)*(2.D0/27.D0)
      CALL F(X+H*(2.E0/27.E0),W12,W2,JIN,LIN,JSIN,MMS)
      DO 201 I=1,N
201   W12(I)=Y(I)+H*(W1(I)*(1.D0/36.D0) + W2(I)*(1.D0/12.D0) )
      CALL F(X+H*(1.D0/9.D0),W12,W3,JIN,LIN,JSIN,MMS)
      DO 202 I=1,N
202   W12(I)=Y(I)+H*(W1(I)*(1.D0/24.D0) + W3(I)*(1.D0/8.D0) )
      CALL F(X+H*(1.D0/6.D0),W12,W4,JIN,LIN,JSIN,MMS)
      DO 203 I=1,N
203   W12(I)=Y(I)+H*(W1(I)*(5.D0/12.D0) - W3(I)*(25.D0/16.D0)
     *                                  + W4(I)*(25.D0/16.D0) )
      CALL F(X+H*(5.D0/12.D0),W12,W5,JIN,LIN,JSIN,MMS)
      DO 204 I=1,N
204   W12(I)=Y(I)+H*(W1(I)*(1.D0/20.D0) + W4(I)*(1.D0/4.D0)
     *                                  + W5(I)*(1.D0/5.D0) )
      CALL F(X+H*0.5D0,W12,W6,JIN,LIN,JSIN,MMS)
      DO 205 I=1,N
205   W12(I)=Y(I)+H*(-W1(I)*(25.D0/108.D0) + W4(I)*(125.D0/108.D0)
     *               -W5(I)*(65.D0/27.D0)  + W6(I)*(125.D0/54.D0) )
      CALL F(X+H*(5.D0/6.D0),W12,W7,JIN,LIN,JSIN,MMS)
      DO 206 I=1,N
206   W12(I)=Y(I)+H*(W1(I)*(31.D0/300.D0) + W5(I)*(61.D0/225.D0)
     *              -W6(I)*(2.D0/9.D0)    + W7(I)*(13.D0/900.D0) )
      CALL F(X+H*(1.D0/6.D0),W12,W8,JIN,LIN,JSIN,MMS)
      DO 207 I=1,N
207   W12(I)=Y(I)+H*(W1(I)*2.D0           - W4(I)*(53.D0/6.D0)
     *             + W5(I)*(704.D0/45.D0) - W6(I)*(107.D0/9.D0)
     *             + W7(I)*(67.D0/90.D0)  + W8(I)*3.D0 )
      CALL F(X+H*(2.D0/3.D0),W12,W9,JIN,LIN,JSIN,MMS)
      DO 208 I=1,N
208   W12(I)=Y(I)+H*(-W1(I)*(91.D0/108.D0)  + W4(I)*(23.D0/108.D0)
     *              - W5(I)*(976.D0/135.D0) + W6(I)*(311.D0/54.D0)
     *              - W7(I)*(19.D0/60.D0)   + W8(I)*(17.D0/6.D0)
     *                                      - W9(I)*(1.D0/12.D0) )
      CALL F(X+H*(1.D0/3.D0),W12,W10,JIN,LIN,JSIN,MMS)
      DO 209 I=1,N
209   W12(I)=Y(I)+H*(W1(I)*(2383.D0/4100.D0) - W4(I)*(341.D0/164.D0)
     *             + W5(I)*(4496.D0/1025.D0) - W6(I)*(301.D0/82.D0)
     *             + W7(I)*(2133.D0/4100.D0) + W8(I)*(45.D0/82.D0)
     *             + W9(I)*(45.D0/164.D0)    + W10(I)*(18.D0/41.D0) )
      CALL F(XTRIAL,W12,W11,JIN,LIN,JSIN,MMS)
      DO 210 I=1,N
210   W12(I)=Y(I)+H*(W1(I)*(3.D0/205.D0) - W6(I)*(6.D0/41.D0)
     *             - W7(I)*(3.D0/205.D0) - W8(I)*(3.D0/41.D0)
     *             + W9(I)*(3.D0/41.D0)  + W10(I)*(6.D0/41.D0) )
      CALL F(X,W12,W2,JIN,LIN,JSIN,MMS)
      DO 211 I=1,N
211   W12(I)=Y(I)+H*(-W1(I)*(1777.D0/4100.D0) - W4(I)*(341.D0/164.D0)
     *              + W5(I)*(4496.D0/1025.D0) - W6(I)*(289.D0/82.D0)
     *              + W7(I)*(2193.D0/4100.D0) + W8(I)*(51.D0/82.D0)
     *              + W9(I)*(33.D0/164.D0)    + W10(I)*(12.D0/41.D0)
     *                                        + W2(I) )
      CALL F(XTRIAL,W12,W3,JIN,LIN,JSIN,MMS)
C
C     ESTIMATE ERROR
C
      EST=0.D0
      DO 300 I=1,N
300   EST=DMAX1(EST,ABS((41.D0/840.D0)*(W1(I)+W11(I)-W2(I)-W3(I))))
C
C     DECISIONS
C
      IF(EST.LT.TOL) GOTO 400
      GOTO 402
400   IND=2
      X=XTRIAL
      DO 401 I=1,N
401   Y(I)=Y(I)+H*(W1(I)*(41.D0/840.D0) + W6(I)*(34.D0/105.D0)
     *           + W7(I)*(9.D0/35.D0)   + W8(I)*(9.D0/35.D0)
     *           + W9(I)*(9.D0/280.D0)  + W10(I)*(9.D0/280.D0)
     *                                  + W11(I)*(41.D0/840.D0) )
      IF(X.EQ.XEND) GOTO 406
      GOTO 405
402   IND=3
      IF(ABS(H).LE.HMIN) GOTO 403
      GOTO 404
403   IND=-2
      HI=0.D0
      RETURN
404   CONTINUE
405   CONTINUE
      GOTO 100
406   HI=SIGN(HMAG,H)
      RETURN
      END




