      SUBROUTINE BVPLSQ (FCN,BC,IVPSOL,NODE,NBC,M,T,X,EPS,ITMAX,
     1                   IFLAG,NONLIN,NRW,RW,NIW,IW,NI2W,I2W)
C*    Begin Prologue BVPLSQ
      EXTERNAL FCN,BC,IVPSOL
      DOUBLE PRECISION EPS,T(M),X(NODE,M),RW(NRW)
      INTEGER IW(NIW)
      INTEGER I2W(NI2W)
C
C---------------------------------------------------------------------
C
C*  Title
C
C     (B)oundary (V)alue (P)roblem (L)east (Sq)uares Solver for highly 
C     nonlinear (possibly overdetermined) two point boundary value
C     problems.
C
C*  Written by        P. Deuflhard, G.Bader
C*  Purpose           Solution of overdetermined nonlinear two-point
C                     boundary value problems.
C*  Method            Local Nonlinear two-point Boundary Value
C                     least squares problems solver
C                     (Multiple shooting approach)
C*  Category          I1b2a - Differential and integral equations
C                             Two point boundary value problems
C*  Keywords          Nonlinear boundary value problems, Multiple
C                     shooting, Gauss Newton methods
C*  Version           0.6 (Test version)
C*  Revision          January 1991
C*  Latest Change     January 1991
C*  Library           CodeLib
C*  Code              Fortran 77, Double Precision
C*  Environment       Standard Fortran 77 environment on PC's,
C                     workstations and hosts.
C*  Copyright     (c) Konrad-Zuse-Zentrum fuer
C                     Informationstechnik Berlin (ZIB)
C                     Takustrasse 7, D-14195 Berlin-Dahlem
C                     phone : + 49/30/84185-0
C                     fax   : + 49/30/84185-125
C*  Contact           Bodo Erdmann
C                     ZIB, Division Scientific Computing, 
C                          Department Numerical Analysis and Modelling
C                     phone : + 49/30/84185-185
C                     fax   : + 49/30/84185-107
C                     e-mail: erdmann@zib.de
C
C*    References:
C
C     /1/ R.Bulirsch:
C         Die Mehrzielmethode zur numerischen Loesung von
C         nichtlinearen Randwertproblemen und Aufgaben der
C         optimalen Steuerung.
C         Carl-Cranz-Gesellschaft: Tech.Rep. (Oct.1971)
C
C     /2/ J.Stoer, R.Bulirsch:
C         Einfuehrung in die Numerische Mathematik II.
C         Berlin, Heidelberg, New York: Springer (1st Ed. 1973)
C
C     /3/ P.Deuflhard:
C         A Modified Newton Method for the Solution of
C         Ill-Conditioned Systems of Nonlinear Equations with
C         Application to Multiple Shooting.
C         Numer. Math. 22, 289-315 (1974)
C
C     /4/ P.Deuflhard:
C         Recent Advances in Multiple Shooting Techniques.
C         (Survey Article including further References)
C         In: I.Gladwell, D.K.Sayers (Ed.): Computational
C         Techniques for Ordinary Differential Equations.
C         Section 10, P.217-272.
C         London, New York: Academic Press (1980)
C
C     /5/ P.Deuflhard, G.Bader:
C         Multiple Shooting Techniques Revisited.
C         Univ. Heidelberg, SFB 123, Tech. Rep. 163 (1982)
C
C     /6/ P.Deuflhard, P.Kunkel:
C         Global Versus Local Linear Equation Solving
C         in Multiple Shooting Techniques
C         Univ. Heidelberg, SFB 123, Tech. Rep. ***
C
C  ---------------------------------------------------------------
C
C* Licence
C    You may use or modify this code for your own non commercial
C    purposes for an unlimited time.
C    In any case you should not deliver this code without a special
C    permission of ZIB.
C    In case you intend to use the code commercially, we oblige you
C    to sign an according licence agreement with ZIB.
C
C* Warranty
C    This code has been tested up to a certain level. Defects and
C    weaknesses, which may be included in the code, do not establish
C    any warranties by ZIB. ZIB does not take over any liabilities
C    which may follow from aquisition or application of this code.
C
C* Software status
C    This code is not under special care of ZIB and belongs to ZIB
C    software class 3.
C
C     ------------------------------------------------------------
C
C  EXTERNAL SUBROUTINES (TO BE SUPPLIED BY THE USER)
C----------------------
C    FCN (T,Z,DZ)            RIGHT-HAND SIDE OF SYSTEM OF FIRST-ORDER
C                            DIFFERENTIAL EQUATIONS
C    BC (XA,XB,R)            TWO-POINT BOUNDARY CONDITIONS
C                            ( A=T(1), B=T(M) )
C    IVPSOL (N,FCN,T,X,TEND,TOL,HMAX,H,KFLAG)  INTEGRATOR
C
C  INPUT PARAMETERS (* MARKS INOUT PARAMETERS)
C------------------
C    NODE          NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS
C    NBC(.GE.NODE) NUMBER OF BOUNDARY CONDITIONS
C    M             NUMBER OF NODES
C                  M.EQ.2    SINGLE SHOOTING
C                  M.GT.2    MULTIPLE SHOOTING
C  * T(M),X(N,M)   INITIAL DATA
C    EPS           REQUIRED RELATIVE PRECISION OF SOLUTION
C    ITMAX         MAXIMUM PERMITTED NUMBER OF ITERATIONS
C  * IFLAG         PRINT PARAMETER
C                  -1        NO PRINT
C                   0        INITIAL DATA
C                            ITERATIVE VALUES OF LEVEL FUNCTIONS
C                            SOLUTION DATA (OR FINAL DATA, RESPECTIVELY)
C                  +1        ADDITIONALLY
C                            ITERATES T(J),X(I,J), I=1,N,J=1,M
C    NONLIN        BVP CLASSIFICATION BY USER
C                   0        LINEAR BVP
C                   1        NONLINEAR BVP
C                            GOOD INITIAL DATA AVAILABLE
C                   2        HIGHLY NONLINEAR BVP
C                            ONLY BAD INITIAL DATA AVAILABLE
C                            SMALL INITIAL DAMPING FACTOR IN
C                            GAUSS NEWTON METHOD
C                   3        HIGHLY NONLINEAR BVP
C                            ONLY BAD INITIAL DATA AVAILABLE
C                            SMALL INITIAL DAMPING FACTOR IN
C                            GAUSS NEWTON METHOD
C                            ADDITIONALLY INITIAL RANK REDUCTION
C                            TO SEPARABLE LINEAR BOUNDARY CONDITIONS
C
C    NRW           DIMENSION OF REAL WORKSPACE RW(NRW)
C                  -----------------------------------------------
C                  NRW .GE. (N*N+1)*M + 10*M*N + NBC*(5*N+8) + 6*N
C                  -----------------------------------------------
C    RW(NRW)       REAL WORKSPACE
C
C    NIW           DIMENSION OF INTEGER WORKSPACE IW(NIW)
C                  ------------------------------------------
C                  NIW .GE. NBC*(2*N+1) + 3*N
C                  ------------------------------------------
C    IW(NIW)       INTEGER WORKSPACE
C
C    NI2W          DIMENSION OF INTEGER WORKSPACE*2 I2W(NI2W)
C                  ------------------------------------------
C                  DUMMY   (FOR COMPATIBILITY WITH BVPSOG)
C                  ------------------------------------------
C    I2W(NI2W)     INTEGER WORKSPACE
C
C  OUTPUT PARAMETERS
C-------------------
C    T(M),X(N,M)   SOLUTION DATA (OR FINAL DATA, RESPECTIVELY)
C    IFLAG         .GT.0     NUMBER OF ITERATIONS PERFORMED
C                            TO OBTAIN THE SOLUTION
C                  .LT.0     BVPLSQ TERMINATION
C                  -1        ITERATION STOPS AT STATIONARY POINT
C                  -2        ITERATION STOPS AFTER ITMAX ITERATION STEPS
C                            (AS INDICATED BY INPUT PARAMETER ITMAX)
C                  -3        INTEGRATOR FAILED
C                            TO COMPLETE THE TRAJECTORY
C                  -4        GAUSS-NEWTON METHOD
C                            FAILED TO CONVERGE
C                  -5        GIVEN INITIAL VALUES INCONSISTENT WITH
C                            SEPARABLE LINEAR BOUNDARY CONDITIONS
C                  -6        ITERATIVE REFINEMENT
C                            FAILED TO CONVERGE
C                  -7        RELIABLE RELATIVE
C                            ACCURACY GREATER THAN 1.D-2
C                  -8        CONDENSING ALGORITHM FOR LINEAR BLOCK
C                            SYSTEM FAILS, USE GLOBAL LINEAR SOLVER
C                            IN BVP ROUTINE BVPSOG
C                  -10       REAL OR INTEGER WORK-SPACE EXHAUSTED
C
C*    End Prologue
      DOUBLE PRECISION DSQRT,EPMACH,FC,RELDIF,SMALL,TOL,XTHR
C
C
      COMMON /MACHIN/ EPMACH, SMALL
      COMMON  /UNIT/  MOUT
      N=NODE
C
C  OUTPUT UNIT FOR ITERATION MONITOR
C------------------------------------
      MOUT=6
C---------------------------------------------------------------------
C       MACHINE DEPENDENT CONSTANTS
C      -----------------------------
C  (ADAPTED TO IBM 3081D, UNIVERSITY OF HEIDELBERG)
C
C  RELATIVE MACHINE PRECISION
      CALL ZIBCONST(EPMACH,SMALL)
C-----------------------------------------------------------------------
C      INTERNAL PARAMETERS
C    -----------------------
C  STANDARD VALUES FIXED BELOW
C
C  SCALING THRESHOLD
      XTHR=SMALL
C
C  PRESCRIBED RELATIVE PRECISION FOR NUMERICAL INTEGRATION
      TOL=EPS*1.D-1
      TOL=EPS*1.D-2
C
C
C  PRESCRIBED RELATIVE DEVIATION FOR NUMERICAL DIFFERENTIATION
      RELDIF=DSQRT(TOL*1.D-1)
C
C  STARTING VALUE OF RELAXATION FACTOR  (1.D-2 .LE. FC .LE. 1.D0)
C  FOR LINEAR OR MILDLY NONLINEAR PROBLEMS
      FC=1.D0
C  FOR HIGHLY NONLINEAR PROBLEMS
      IF(NONLIN.GT.1) FC=1.D-2
C
C  STARTING VALUE FOR PSEUDO-RANK OF SENSITIVITY MATRIX E
      IRANK=N
C
C---------------------------------------------------------------------
C
C  CHECK FOR SUFFICIENT REAL/INTEGER WORKSPACE
C----------------------------------------------
      MINRW=(N*N+1)*M + 10*M*N + NBC*(5*N+8) + 6*N
      MINIW=NBC*(2*N+1) + 3*N
      IF(IFLAG.GE.0) WRITE(MOUT,1000) MINRW,MINIW
      IF(MINRW.GT.NRW .OR. MINIW.GT.NIW) GOTO 900
C
C  INITIAL PREPARATIONS
C
      M1=M-1
      NN=N*N
      NM=N*M
      NBCN=NBC*N
      NM1=N*M1
C
      I2=NBC+1
      I3=I2+N
      I4=I3+N
      I5=I4+N
      I6=I5+NBCN
C
      N2=NN*M1+1
      N3=N2+NBCN
      N4=N3+NBCN
      N5=N4+NBCN
      N6=N5+NBCN
      N7=N6+NN
      N8=N7+NM
      N9=N8+NM
      N10=N9+NM
      N11=N10+NM
      N12=N11+NM
      N13=N12+NM
      N14=N13+NM1
      N15=N14+NM1
      N16=N15+NM1
      N17=N16+NM1
      N18=N17+N
      N19=N18+NBC
      N20=N19+NBC
      N21=N20+NBC
      N22=N21+NBC
      N23=N22+NBC
      N24=N23+NBC
      N25=N24+N
      N26=N25+N
      N27=N26+N
      N28=N27+N
      N29=N28+NBC
      N30=N29+N
      N31=N30+M
      N32=N31+NBC
C
      CALL BFBVPL(FCN,BC,IVPSOL,N,NBC,M,M1,T,X,EPS,TOL,RELDIF,FC,NONLIN
     1  ,IRANK,ITMAX,IFLAG, IW(1),IW(I2),IW(I3),IW(I4),IW(I5),IW(I6)
     2  ,RW(1),RW(N2),RW(N3),RW(N4),RW(N5),RW(N6),RW(N7),RW(N8),RW(N9)
     3  ,RW(N10),RW(N11),RW(N12),RW(N13),RW(N14),RW(N15),RW(N16),RW(N17)
     4  ,RW(N18),RW(N19),RW(N20),RW(N21),RW(N22),RW(N23),RW(N24),RW(N25)
     5  ,RW(N26),RW(N27),RW(N28),RW(N29),RW(N30),RW(N31),RW(N32),XTHR)
C
C  SOLUTION EXIT
      RETURN
C
C  FAIL EXIT  WORK-SPACE EXHAUSTED
900   IF(IFLAG.GE.0.AND.MINRW.GT.NRW) WRITE(MOUT,1001)
      IF(IFLAG.GE.0.AND.MINIW.GT.NIW) WRITE(MOUT,1002)
      IFLAG=-10
      RETURN
C
1000  FORMAT(30H0 MINIMAL REQUIRED WORK-SPACE:,/,
     @       20H0 REAL    ARRAY  RW(,I4,1H)/
     @       20H0 INTEGER ARRAY  IW(,I4,1H))
1001  FORMAT(35H0 ERROR:  REAL WORK-SPACE EXHAUSTED,/)
1002  FORMAT(38H0 ERROR:  INTEGER WORK-SPACE EXHAUSTED,/)
C
C  END DRIVER ROUTINE BVPLSQ
C
      END
C
      SUBROUTINE BFBVPL (FCN,BC,IVPSOL,N,NBC,M,M1,T,X,EPS,TOL,RELDIF,FC
     1           ,NONLIN,IRANK,ITMAX,IFLAG, IROW,ICOL,ICOLB,PIVOT,IA,IB
     2           ,G,A,B,BG,E,QE,DX,DDX,DXQ,DXQA,XA,XW,XU,HH,DHH,HHA
     3           ,D,DE,R,DR,RA,U,DU,QU,X1,XM,T1,T2,DX1,RF,US,EH,XTHR)
C
C
C---------------------------------------------------------------------
C         SOLUTION OF (POSSIBLY) OVERDETERMINED NONLINEAR
C                TWO-POINT BOUNDARY VALUE PROBLEMS
C---------------------------------------------------------------------
C
       DOUBLE PRECISION T(M),X(N,M)
     1 ,G(N,N,M1), A(NBC,N), B(NBC,N), BG(NBC,N), E(NBC,N), QE(N,N)
     2 ,DX(N,M), DDX(N,M), DXQ(N,M),  DXQA(N,M), XA(N,M), XW(N,M)
     3 ,XU(N,M1), HH(N,M1), DHH(N,M1), HHA(N,M1)
     4 ,D(N), DE(NBC),  R(NBC), DR(NBC), RA(NBC), U(NBC), DU(NBC), QU(N)
     5 ,X1(N), XM(N), T1(N), T2(NBC), DX1(N), RF(M)
     6 ,US(NBC), EH(NBC,N)
C
      DOUBLE PRECISION     COND  ,COND1 ,COND2 , CORR ,CONV  ,CONVA ,
     1DABS  ,DEL   ,DSQRT ,EPH   ,EPMACH,EPS   ,EPSMIN,FC    ,
     2FCA   ,FCDNM ,FCH   ,FCMIN ,FCMINH,FCMIN2,FCNUM ,FCS   ,EIGHT ,
     3H     ,HALF  ,HFIN  ,HMAX  ,HSAVE ,HSTART,ONE   ,REDH  ,RELDIF,
     4S     ,SENS1 ,SENS2 ,SIGDEL,SIGDLH,SIGMA ,SKAP  ,SMALL ,ST    ,
     5SUMF  ,SUMX  ,SUMXA ,TEN   ,TFAIL ,TH    ,TJ    ,TJ1   ,TOL   ,
     6TOLH  ,TOLMIN,TWO   ,ZERO  ,EPX1  ,EPDX1 ,EPX1H ,XTHR  ,CONDE ,
     7DMAX1
C
      INTEGER IROW(NBC),ICOL(N),ICOLB(N),PIVOT(N),IA(NBC,N),IB(NBC,N)
C
      COMMON /MACHIN/ EPMACH, SMALL
      COMMON  /UNIT/  MOUT
C
      EXTERNAL FCN,BC,IVPSOL,BFSCLE
C
      DATA  REDH/1.D-2/,ZERO/0.D0/,HALF/0.5D0/,
     @FCS/0.7D0/ , ONE/1.D0/ , TWO/2.D0/ , EIGHT/8.D0/ , TEN/1.D1/
C
C---------------------------------------------------------------------
C      INTERNAL PARAMETERS
C    -----------------------
C  STANDARD VALUES FIXED BELOW
C
C  MAXIMUM PERMITTED SUB-CONDITION NUMBER OF MATRIX E
      CONDE=ONE/EPMACH
C
C  PRESCRIBED RELATIVE PRECISION IN BVPLSQ
      EPSMIN=DSQRT(TEN*EPMACH)
      IF(EPS.LT.EPSMIN) EPS=EPSMIN
C
C  MINIMUM RELATIVE PRECISION OF INTEGRATOR (TO BE ADAPTED)
      TOLMIN=EPMACH*TEN*TEN
      IF(TOL.LE.ZERO) TOL=EPS/TEN
      IF(TOL.LT.TOLMIN) TOL=TOLMIN
C
C
C  MINIMUM PERMITTED VALUE OF RELAXATION FACTOR
      FCMIN=1.D-2
      IF(FC.LT.FCMIN) FC=FCMIN
      IF(FC.GT.ONE) FC=ONE
C
C
C  MAXIMUM PERMITTED NUMBER OF ITERATIVE REFINEMENTS SWEEPS
      NYMAX=M-1
      IF(NBC.GT.N) NYMAX=0
C
C
C  DECISION PARAMETER FOR JACOBIAN RANK-1 UPDATES (SIGMA.GT.1.)
C  RANK-1 UPDATES INHIBITED, IF SIGMA.GT.1./FCMIN IS SET
      SIGMA=2.D0
C
C---------------------------------------------------------------------
C
C
C  INITIAL PREPARATIONS
C----------------------
      KPRINT=IFLAG
      FCMIN2=FCMIN*FCMIN
      FCMINH=DSQRT(FCMIN)
      ITER=0
      KOUNT=0
      INIT=0
      LEVEL=0
      IREPET=0
      INUM=1
      IRKMAX=0
      FCA=FC
      HSTART=(T(2)-T(1))*REDH
      SENS1=ZERO
      SENS2=ZERO
      COND1=ONE
      COND2=ONE
      SIGDLH=ZERO
      IRANKB=0
      SUMF=ONE
C
      DO 100 I=1,NBC
      DO 100 K=1,N
      IA(I,K)=0
100   IB(I,K)=0
      DO 101 I=1,N
      DO 101 K=1,N
101   G(I,K,1)=ZERO
C
      IF(KPRINT.LT.0) GOTO 2000
      WRITE(MOUT,1001)
      WRITE(MOUT,1002)
      DO 200 J=1,M
200   WRITE(MOUT,1003) T(J),(X(I,J),I=1,N)
      WRITE(MOUT,1004) N,NBC,M,EPS,ITMAX
      WRITE(MOUT,1001)
      IF(KPRINT.GT.0) GOTO 2000
      WRITE(MOUT,1005)
      WRITE(MOUT,1006)
      GOTO 2000
C
C
C  PRELIMINARY NEW ITERATE
C-------------------------
1000  INIT=1
      DO 1100 J=1,M
      DO 1100 I=1,N
1100  X(I,J)=XA(I,J)+FC*DX(I,J)
      IF(ITER.GE.ITMAX) GOTO 9200
C
C
C  COMPUTATION OF THE TRAJECTORIES
C  (SOLUTION OF M1 INITIAL VALUE PROBLEMS)
C-----------------------------------------
2000  J=1
      KOUNT=KOUNT+1
      HSAVE=HSTART
2100  J1=J+1
      TJ=T(J)
      TJ1=T(J1)
      H=HSAVE
      HMAX=DABS(TJ1-TJ)
      KFLAG=0
      DO 2110 K=1,N
2110  T1(K)=X(K,J)
      CALL IVPSOL (N,FCN,TJ,T1,TJ1,TOL,HMAX,H,KFLAG)
      HFIN=HSAVE
      HSAVE=H
      IF(H.NE.ZERO) GOTO 2200
C
C  SINGULAR TRAJECTORY
      IFAIL=TJ
      KFLAG=-J
      IF(INIT.EQ.0) GOTO 9300
      IF(KPRINT.GE.0) WRITE(MOUT,2001)
      FC=FC*HALF
      IF(FC.LT.FCMIN) GOTO 7700
      GOTO 1000
C
C  CONTINUITY CONDITIONS
C-----------------------
2200  DO 2210 K=1,N
      TH=T1(K)
      XU(K,J)=TH
2210  HH(K,J)=TH-X(K,J1)
      J=J1
      IF(J.LT.M) GOTO 2100
C
C  TWO-POINT BOUNDARY CONDITIONS
C-------------------------------
      DO 2300 I=1,N
      XM(I)=X(I,M)
2300  X1(I)=X(I,1)
      CALL BC (X1,XM,R)
C
      IF(INIT.EQ.0) GOTO 5100
      LEVEL=1
C
C  COMPUTATION OF CONDENSED RIGHT-HAND SIDE U(NRE)
C------------------------------------------------
3000  IF(IRANK.GT.0) CALL BFRHS1 (N,NBC,NRE,NE,M1,1,HH,R,B,G,U,DE,
     @                          T1,BG,IROW)
C
C  SAVING OF RIGHT-HAND SIDE U
      IF (IRANK.LT.NE) GOTO 3020
      DO 3010 I=1,NRE
3010  US(I)=U(I)
3020  CONTINUE
C
C  (BEST) CONSTRAINED LEAST SQUARES SOLUTION OF LINEAR (NRE,NE)-SYSTEM
C--------------------------------------------------------------------
3100  DO 3110 I=1,N
3110  DX1(I)=ZERO
C
      IF(IRANK.GT.0)  CALL BVSOLC
     @      (E,NBC,N,IRANKB,NRE,NE,DX1,U,IRANK,D,PIVOT,IREPET,QE,T1)
C
      IF(LEVEL.GT.0 .OR. IREPET.NE.0 .OR. IRANK.EQ.0) GOTO 3130
      DO 3120 I=1,IRANK
3120  QU(I)=U(I)
3130  CONTINUE
C
C  ITERATIVE REFINEMENT OF DX1
      EPH=EPS
      IF(IRANK.LT.NE .OR. NE.EQ.0) GOTO 3160
      DO 3140 I=1,NRE
      S=ZERO
      DO 3141 K=1,NE
3141  S=S+EH(I,K)*DX1(K)
3140  DU(I)=US(I)-S
C  SOLUTION OF RESIDUAL EQUATION
      CALL BVSOLC (E,NBC,N,IRANKB,NRE,NE,T2,DU,IRANK,D,PIVOT,
     @             IREPET,QE,T1)
      EPX1=ZERO
      EPDX1=ZERO
      DO 3150 I=1,NE
      S=DABS(T2(I))
      IF(EPDX1.LT.S) EPDX1=S
      DX1(I)=DX1(I)+T2(I)
      S=DABS(DX1(I))
      IF(EPX1.LT.S) EPX1=S
3150  CONTINUE
      EPX1H=EPDX1/EPX1
      EPH=TEN*EPDX1
      IF(EPX1H.GT.HALF) GOTO 9800
3160  CONTINUE
C
C  DESCALING AND BACK-PERMUTATION OF SOLUTION DX1
      DO 3170 L=1,N
      I=ICOL(L)
3170  DXQ(I,1)=DX1(L)*XW(I,1)
C
C
C  RECURSIVE COMPUTATION OF DXQ(N,2),...,DXQ(N,M)
C-------------------------------------------------
      CALL BFRECU (N,M,M1,1,HH,G,DXQ,T1,T2)
C
C
C         ITERATIVE REFINEMENT SWEEPS  NY=1,..,NYMAX
C-----------------------------------------------------------------------
      NY=0
      SIGDEL=TEN
      SIGDEL=TEN*TEN
      SIGDLH=ZERO
      IF(EPH.LT.EPS) EPH=EPS
      IF(NYMAX.EQ.0 .OR. IRANK.LT.NE .OR. NE.EQ.0) GOTO 3990
      IF(KPRINT.GT.0) WRITE(MOUT,3201)
C
C  COMPUTATION OF REQUIRED CONTINUITY RESIDUALS DHH(N,M1)
      JN=1
      JIN=M
      GOTO 3230
C
3200  DO 3220 J=JN,M1
      J1=J+1
      DO 3220 I=1,N
      S=HH(I,J)
      DO 3221 K=1,N
3221  S=S+G(I,K,J)*DXQ(K,J)
3220  DHH(I,J)=S-DXQ(I,J1)
C
C  COMPUTATION OF BOUNDARY RESIDUAL DR(NBC)
3230  DO 3240 I=1,NBC
      S=R(I)
      DO 3241 K=1,N
3241  S=S+A(I,K)*DXQ(K,1)+B(I,K)*DXQ(K,M)
3240  DR(I)=S
C
C  COMPUTATION OF CONDENSED RESIDUAL DU(NRE)
      IF(IRANK.GT.0) CALL BFRHS1 (N,NBC,NRE,NE,M1,JIN,DHH,DR,B,G,DU,DE,
     @                          T1,BG,IROW)
C
C  COMPUTATION OF CORRECTION DDX(N,1)
      DO 3250 I=1,N
3250  DX1(I)=ZERO
      IF(IRANK.GT.0) CALL BVSOLC
     @       (E,NBC,N,IRANKB,NRE,NE,DX1,DU,IRANK,D,PIVOT,IREPET,QE,T1)
C
C  DESCALING OF DDX(N,1), REFINEMENT OF DXQ(N,1)
C
      CORR=ZERO
      DO 3260 L=1,N
      I=ICOL(L)
      S=DX1(L)
      IF(CORR.LT.DABS(S)) CORR=DABS(S)
      S=S*XW(I,1)
      DDX(I,1)=S
3260  DXQ(I,1)=DXQ(I,1)+S
      IF (CORR.LT.EPH) GOTO 3261
      EPH=CORR
      GOTO 9800
3261  RF(1)=CORR
C
C  RECURSIVE COMPUTATION OF DDX(N,2),...,DDX(N,M)
      CALL BFRECU (N,M,M1,JIN,DHH,G,DDX,T1,T2)
C
C  REFINEMENT OF DXQ(N,2),...,DXQ(N,M)
      DO 3270 J=2,M
      CORR=ZERO
      DO 3271 I=1,N
      S=DDX(I,J)
      DXQ(I,J)=DXQ(I,J)+S
      S=DABS(S)/XW(I,J)
      IF(CORR.LT.S) CORR=S
3271  CONTINUE
3270  RF(J)=CORR
C
C  DETERMINATION OF SWEEP INDEX JN
      JA=JN
      DO 3280 J=1,M
      IF(RF(J).GT.EPH) GOTO 3290
3280  JN=J
C
3290  NY=NY+1
      IF(KPRINT.GT.0) WRITE(MOUT,3202) NY,JA,(RF(J),J=1,M)
      IF(JN.LE.JA) GOTO 9600
      IF(JN.EQ.M) GOTO 3900
      JIN=JN
      IF(NY.GT.1 .OR. LEVEL.EQ.0) GOTO 3200
C
C  DETERMINATION AND ADAPTATION OF PARAMETERS TOL AND RELDIF
3900  IF(LEVEL.EQ.0 .OR. NY.GT.1) GOTO 3920
      DO 3910 J=1,M1
      S=RF(J+1)/RF(J)
      IF(SIGDLH.LT.S) SIGDLH=S
      RF(J)=S
3910  CONTINUE
      IF(KPRINT.GT.0) WRITE(MOUT,3205) (RF(J),J=1,M1)
      SIGDEL=DMAX1(SIGDLH,SIGDEL)
      TH=TOL*SIGDEL
C
      IF(TH.GT.REDH) GOTO 9700
C
      IF(TH.GT.EPH) EPH=TH
      TOLH=EPS/SIGDEL
      IF(TOLH.LT.TOLMIN) TOLH=TOLMIN
      TOL=TOLH
      RELDIF=DSQRT(TOL/SIGDEL)
      IF(KPRINT.LT.0) GOTO 3920
      WRITE(MOUT,3203) TOLH,RELDIF
      WRITE(MOUT,3204)
3920  IF(JN.NE.M) GOTO 3200
C
C      END OF ITERATIVE REFINEMENT SWEEPS
C-----------------------------------------------------------------------
C
C
C
C  PROJECTION OF SEPARATED LINEAR BOUNDARY CONDITIONS AT T(M)
3990  IF(NB.EQ.0) GOTO 4000
      DO 3991 K=1,NB
      IC=ICOLB(K)
3991  DXQ(IC,M)=ZERO
C
C  EVALUATION OF SCALED NATURAL LEVEL FUNCTION SUMX
C  AND SCALED MAXIMUM ERROR NORM CONV
C---------------------------------------------------
4000  CONV=ZERO
      SUMX=ZERO
      DO 4010 J=1,M
      DO 4010 I=1,N
      S=DABS(DXQ(I,J))/XW(I,J)
      IF(CONV.LT.S) CONV=S
4010  SUMX=SUMX+S*S
C
C  EVALUATION OF SCALED STANDARD LEVEL FUNCTION SUMF
C---------------------------------------------------
      SUMF=ZERO
      DO 4020 I=1,NBC
4020  SUMF=SUMF + (R(I)*DE(I)/SMALL)**2
      DO 4030 J=1,M1
      J1=J+1
      DO 4030 I=1,N
4030  SUMF=SUMF + (HH(I,J)/XW(I,J1))**2
C
      IF(LEVEL.GT.0) GOTO 4500
C
C  ORDINARY GAUSS-NEWTON CORRECTIONS DX(N,M)
C--------------------------------------------
      SUMXA=SUMX
      CONVA=CONV
      DO 4110 J=1,M
      DO 4110 I=1,N
      XA(I,J)=X(I,J)
4110  DX(I,J)=DXQ(I,J)
C
C  A-PRIORI ESTIMATE OF RELAXATION FACTOR FC
C--------------------------------------------
4200  JRED=0
      IF(ITER.EQ.0 .OR. NONLIN.EQ.0) GOTO 4400
      IF( (NEW.GT.0 .OR. IRANK.LT.NE.AND.IRANKA.LT.NE)
     @                              .AND.IREPET.EQ.0 ) GOTO 4350
C
C  FULL RANK CASE (INDEPENDENT OF PRECEDING RANK)
C  COMPUTATION OF THE DENOMINATOR OF A-PRIORI ESTIMATE
      FCDNM=ZERO
      DO 4201 J=1,M
      DO 4201 I=1,N
4201  FCDNM=FCDNM+((DX(I,J)-DXQA(I,J))/XW(I,J))**2
C
C  RANK-DEFICIENT CASE (IF PREVIOUS RANK WAS FULL)
C  COMPUTATION OF THE PROJECTED DENOMINATOR OF A-PRIORI ESTIMATE
      IF(IRANK.EQ.NE) GOTO 4300
      DO 4210 L=1,N
4210  T1(L)=ZERO
      DEL=ZERO
      IF(IRANK.EQ.0) GOTO 4240
C
      DO 4220 L=1,NE
      K=ICOL(L)
4220  DX1(L)=DXQA(K,1)/XW(K,1)
C
C  PROJECTION FOR REDUCED COMPONENT DX1(NE)
      CALL BFPRJT (N,NE,IRANK,DEL,DX1,D,T2,QE,PIVOT)
C
      DO 4230 L=1,NE
      K=ICOL(L)
4230  T1(K)=DX1(L)*XW(K,1)
C
4240  DO 4260 J=1,M1
      DO 4250 I=1,N
      S=ZERO
      DO 4251 K=1,N
4251  S=S+T1(K)*G(I,K,J)
      J1=J+1
      ST=ONE/XW(I,J1)
      DEL=DEL+S*ST*ST*(S+(DX(I,J1)-DXQA(I,J1))*TWO)
4250  T2(I)=S
      DO 4260 I=1,N
4260  T1(I)=T2(I)
      FCDNM=FCDNM+DEL
C
4300  FC=FCA/FCMIN
      IF(FCDNM.GT.FCNUM*FCMIN2) FC=DSQRT(FCNUM/FCDNM)*FCA
C
4350  IREPET=0
      IF(FC.LT.FCMIN) GOTO 7700
      IF(FC.GT.FCS) FC=ONE
C
4400  IF(KPRINT.LT.0) GOTO 1000
      IF(KPRINT.GT.0) WRITE(MOUT,1005)
      IF(KPRINT.GT.0) WRITE(MOUT,1006)
      IRKDEF=NE-IRANK
      WRITE(MOUT,4401) ITER,NY,SUMF,SUMXA,NEW,IRKDEF
      IF(KPRINT.GT.0) WRITE(MOUT,1005)
      GOTO 1000
C
C  SIMPLIFIED GAUSS-NEWTON CORRECTIONS DXQ(N,M)
C-----------------------------------------------
C
C  RANK INDEPENDENT CONVERGENCE TEST
C
4500  IF(CONV.LE.EPS.AND.IRKMAX.EQ.NE) GOTO 9000
C
C  NATURAL MONOTONICITY TEST
C
      IF(SUMX.LE.SUMXA) GOTO 5000
C
C  REDUCTION OF RELAXATION FACTOR FC
C------------------------------------
      IF(KPRINT.LT.0) GOTO 4600
      IF(KPRINT.EQ.0) GOTO 4610
      WRITE(MOUT,1005)
      WRITE(MOUT,1006)
4610  WRITE(MOUT,5001) ITER,NY,SUMF,SUMX,FC
      IF(KPRINT.GT.0) WRITE(MOUT,1005)
C
4600  JRED=JRED+1
      IF(NONLIN.EQ.0) GOTO 9800
      TH=DSQRT(SUMX/SUMXA)
      TH=DSQRT(EIGHT*(TH+FC-ONE)/FC+ONE)-ONE
      FC=FC/TH
      IF(FC.LT.FCMIN .OR. NEW.GT.0.AND.JRED.GT.1) GOTO 7700
      GOTO 1000
C
C  PREPARATIONS TO START THE FOLLOWING ITERATION STEP
C----------------------------------------------------
5000  ITER=ITER+1
      LEVEL=0
C
      IF(KPRINT.LT.0) GOTO 5100
      IF(KPRINT.EQ.0) GOTO 5010
      WRITE(MOUT,1005)
      WRITE(MOUT,1006)
5010  WRITE(MOUT,5001) ITER,NY,SUMF,SUMX,FC
      IF(KPRINT.EQ.0) GOTO 5100
      WRITE(MOUT,1005)
      DO 5020 J=1,M
5020  WRITE(MOUT,1003) T(J),(X(I,J),I=1,N)
C
C
C  SCALING OF VARIABLES X(N,M)
5100  CALL BFSCLE (N,M,M1,X,XU,XW,XTHR)
C
      IF(INIT.EQ.0) GOTO 6000
C
C  SAVING OF VALUES DXQ(N,M)
      DO 5200 J=1,M
      DO 5200 I=1,N
5200  DXQA(I,J)=DXQ(I,J)
C
C  PRELIMINARY PSEUDO-RANK
      IRANKA=IRANK
      IF(IRANK.GE.NB.AND.FC.GT.FCMINH) IRANK=NE
      IF(IRANK.LT.NB) IRANK=NB
C
C  A-POSTERIORI ESTIMATE OF RELAXATION FACTOR FC
C------------------------------------------------
      FCNUM=ZERO
      FCDNM=ZERO
      TH=FC-ONE
      DO 5400 J=1,M
      DO 5400 I=1,N
      FCNUM=FCNUM+(DX(I,J)/XW(I,J))**2
5400  FCDNM=FCDNM+((DXQ(I,J)+TH*DX(I,J))/XW(I,J))**2
      FCH=DSQRT(FCNUM/FCDNM)*FC*FC*HALF
C
C  DECISION CRITERION FOR JACOBIAN UPDATING TECHNIQUE
C  INUM=0: RANK1 UPDATING, INUM=1: NUMERICAL DIFFERENTIATION
      INUM=0
      IF(FC.LT.FCA.AND.NEW.GT.0 .OR. FCH.LT.FC*SIGMA
     @     .OR. EPH*REDH.GT.EPS .OR. IRANK.GT.IRANKA)  INUM=1
      FCA=FC
      IF(NONLIN.GT.0) FC=FCH
C
C
C DIFFERENCE APPROXIMATION OF BOUNDARY DERIVATIVE MATRICES A AND B
C------------------------------------------------------------------
6000  CONTINUE
      CALL BFDERA (BC,N,NBC,M,XW,X1,XM,R,T2,A,B,RELDIF)
C
C  DETERMINATION OF SPARSE STRUCTURE OF MATRICES A AND B
C  AND DETERMINATION OF INTERNAL ROW SCALING OF SENSITIVITY MATRIX E
C--------------------------------------------------------------------
      ISUM=0
      DO 6100 I=1,NBC
      S=ZERO
      DO 6110 K=1,N
      TH=DABS(A(I,K))*XW(K,1)
      IF(S.LT.TH) S=TH
      TH=DABS(B(I,K))*XW(K,M)
      IF(S.LT.TH) S=TH
      IF(IA(I,K).GT.0) GOTO 6111
      IF(A(I,K).EQ.ZERO) GOTO 6111
      IA(I,K)=1
      ISUM=1
6111  IF(IB(I,K).GT.0) GOTO 6110
      IF(B(I,K).EQ.ZERO) GOTO 6110
      IB(I,K)=1
      ISUM=1
6110  CONTINUE
      IF(S.LT.XTHR) S=XTHR
6100  DE(I)=SMALL/S
C
      IF(ISUM.EQ.0) GOTO 6400
C
C  DETERMINATION OF ROW AND COLUMN PERMUTATION VECTORS
      DO 6210 I=1,N
      ICOL(I)=I
6210  ICOLB(I)=I
      DO 6211 I=1,NBC
6211  IROW(I)=I
C
C  SEARCH FOR SEPARABLE LINEAR BOUNDARY CONDITIONS AT T(1)
      NE=N
      NRE=NBC
      DO 6220 I=1,NBC
      DO 6221 K=1,N
      IF(IB(I,K).NE.0) GOTO 6220
6221  CONTINUE
      ISUM=0
      DO 6222 K=1,N
      IF(IA(I,K).EQ.0) GOTO 6222
      ISUM=ISUM+1
      ICA=K
6222  CONTINUE
      IF(ISUM.GT.1) GOTO 6220
      DO 6223 IS=1,N
      IH=ICOL(IS)
      IF(IH.EQ.ICA) ICOL(IS)=ICOL(NE)
6223  CONTINUE
      DO 6224 IS=1,NBC
      IH=IROW(IS)
      IF(IH.EQ.I) IROW(IS)=IROW(NE)
6224  CONTINUE
      ICOL(NE)=ICA
      NE=NE-1
      IROW(NRE)=I
      NRE=NRE-1
      IF(DABS(R(I)).GT.TEN*EPMACH*DABS(X(ICA,1))) GOTO 9500
6220  CONTINUE
      IF(KPRINT.GE.0.AND.NE.EQ.0) WRITE(MOUT,6001)
      IF(IRANK.GT.NE) IRANK=NE
      IRANKA=IRANK
C
C  SEARCH FOR SEPARABLE LINEAR BOUNDARY CONDITIONS AT T(M)
      NB=0
      IF(NE.EQ.0) GOTO 6400
      DO 6230 I=1,NRE
      IR=IROW(I)
      DO 6231 K=1,N
      IF(IA(IR,K).NE.0) GOTO 6230
6231  CONTINUE
      ISUM=0
      DO 6232 K=1,N
      IF(IB(IR,K).EQ.0) GOTO 6232
      ISUM=ISUM+1
      ICB=K
6232  CONTINUE
      IF(ISUM.GT.1) GOTO 6230
      NB=NB+1
      DO 6233 IS=1,N
      IH=ICOLB(IS)
      IF(IH.EQ.ICB) ICOLB(IS)=ICOLB(NB)
6233  CONTINUE
      ICOLB(NB)=ICB
      IROW(I)=IROW(NB)
      IROW(NB)=IR
      IF(DABS(R(IR)).GT.TEN*EPMACH*DABS(X(ICB,M))) GOTO 9500
6230  CONTINUE
      IF(KPRINT.GE.0.AND.NB.EQ.N) WRITE(MOUT,6001)
C
C  INITIAL RANK STRATEGY FOR HIGHLY NONLINEAR PROBLEMS
      IF(NB.LT.NE.AND.INIT.EQ.0.AND.NONLIN.GT.2) IRANK=NB
C
6400  IF(INUM.EQ.0) GOTO 6800
C
C  DIFFERENCE APPROXIMATION OF WRONSKIAN MATRICES G(1),...,G(M1)
C----------------------------------------------------------------
6500  NEW=0
      KFLAG=0
      CALL BFDERG (FCN,N,NE,M,M1,T,X,XU,XW,T1,TFAIL,G,
     @               ICOL,IVPSOL,HSTART,TOL,RELDIF,KFLAG)
      IF(KFLAG.LT.0) GOTO 9310
C
      IF(M.GT.2) KOUNT=KOUNT+N
      IF(M.EQ.2) KOUNT=KOUNT+NE
C
      GOTO 7000
C
C  RANK-1 UPDATES OF WRONSKIAN MATRICES G(1),...,G(M1)
C------------------------------------------------------
6800  NEW=NEW+1
      CALL BFRK1G (N,M,M1,XW,DX,HH,HHA,T1,G,FCA)
C
C
C
C  COMPUTATION OF (NRE,NE)-SENSITIVITY MATRIX E=-(A+B*G(M1)*...*G(1))
C  (PROJECTIONS INCLUDED)
C-----------------------------------------------------------
C
7000  IF(IRANK.EQ.0) GOTO 7500
      DO 7100 I=1,NRE
      IR=IROW(I)
      DO 7100 K=1,N
7100  E(I,K)=B(IR,K)*DE(IR)
      DO 7200 JJ=1,M1
      J=M-JJ
      DO 7200 I=1,NRE
      DO 7210 K=1,N
      S=ZERO
      DO 7211 L=1,N
7211  S=S+E(I,L)*G(L,K,J)
7210  T1(K)=S
      DO 7220 K=1,N
7220  E(I,K)=T1(K)
7200  CONTINUE
C  INTERNAL ROW AND COLUMN SCALING AND PERMUTATION OF MATRIX E
      DO 7400 K=1,NE
      KC=ICOL(K)
      S=XW(KC,1)
      DO 7400 I=1,NRE
      IR=IROW(I)
7400  E(I,K)=-(A(IR,KC)*DE(IR)+E(I,KC))*S
C
C  SAVE MATRIX E ON EH
      DO 7410 K=1,NE
      DO 7410 I=1,NRE
7410  EH(I,K)=E(I,K)
C
C  RESTORE RANK OF SEPARABLE LINEAR BOUNDARY CONDITIONS AT T(M)
7500  IRANKB=0
      IF(NB.LT.NE) IRANKB=NB
C
C  MONITOR FOR ACTUALLY APPLIED MAXIMUM RANK
      IF(IRKMAX.LT.IRANK) IRKMAX=IRANK
C
C  SAVE VALUES OF  R(N),HH(N,M1)
      IF(IREPET.NE.0) GOTO 7600
      DO 7501 I=1,NBC
7501  RA(I)=R(I)
      DO 7502 I=1,N
      DO 7502 J=1,M1
7502  HHA(I,J)=HH(I,J)
C
C  CONSTRAINED QR-DECOMPOSITION OF (NRE,NE)-MATRIX E
C---------------------------------------------------
7600  COND=CONDE
      IF(IRANK.GT.0) CALL BVDECC
     @            (E,NBC,N,IRANKB,NRE,NE,IRANK,COND,D,PIVOT,
     @            IREPET,QE,T1)
      IF (NONLIN.EQ.0 .AND. IRANK.LT.NE) GOTO 9800
C
C
C  EVALUATION OF SUBCONDITION AND SENSITIVITY NUMBERS
      COND1=ONE
      COND2=ONE
      SENS1=ZERO
      SENS2=ZERO
      IF(IRANKB.EQ.0) GOTO 7650
      SENS1=DABS(D(1))
      COND1=SENS1/DABS(D(IRANKB))
7650  IF(IRANKB.EQ.IRANK) GOTO 7690
      SENS2=DABS(D(IRANKB+1))
      COND2=SENS2/DABS(D(IRANK))
7690  CONTINUE
C
      IF(FCA.LT.1.D0 .OR. FC.LT.1.D0 .OR. ITER.EQ.0) GOTO 7699
      IF(IRANKB.NE.IRANK .AND. SENS2.LT.(EPS/REDH)*SMALL) IRANK=IRANKB
      IF(IRANKB.NE.0 .AND. SENS1.LT.(EPS/REDH)*SMALL) IRANK=0
7699  IF(IREPET) 3100,3000,3000
C
C  RESTORE FORMER VALUES
C------------------------
7700  IREPET=1
      LEVEL=0
      DO 7710 I=1,N
      X1(I)=XA(I,1)
      XM(I)=XA(I,M)
      X(I,1)=XA(I,1)
      DO 7710 J=1,M1
      J1=J+1
      X(I,J1)=XA(I,J1)
      XU(I,J)=X(I,J1)+HHA(I,J)
7710  HH(I,J)=HHA(I,J)
      DO 7711 I=1,NBC
7711  R(I)=RA(I)
      IRKDEF=NE-IRANK
      IF(KPRINT.GE.0) WRITE(MOUT,7701) ITER,FC,IRKDEF
      IF(ITER.EQ.0) FC=FCMIN
      IF(NEW.GT.0) GOTO 6500
C
C  PSEUDO-RANK REDUCTION
C-----------------------
      IREPET=-1
      IF(IRANK.EQ.0) GOTO 9400
      DO 7720 I=1,IRANK
7720  U(I)=QU(I)
      IRANK=IRANK-1
      IF(IRANKB.GT.IRANK) IRANKB=IRANK
      GOTO 7600
C
C
C
C------------------------ EXIT ---------------------------------------
C
C  SOLUTION EXIT
C---------------
C
9000  ITER=ITER+1
      DO 901 J=1,M
      DO 901 I=1,N
901   X(I,J)=X(I,J)+DXQ(I,J)
C
      IF(IRANK.LT.NE.AND.KPRINT.LT.0) GOTO 9100
      IFLAG=ITER
      IF(KPRINT.LT.0) RETURN
      IF(KPRINT.EQ.0) GOTO 9010
      WRITE(MOUT,1005)
      WRITE(MOUT,1006)
9010  WRITE(MOUT,5001) ITER,NY,SUMF,SUMX,FC
      WRITE(MOUT,1005)
      WRITE(MOUT,1001)
      IF(IRANK.LT.NE) GOTO 9100
      WRITE(MOUT,9001) ITER,KOUNT
      WRITE(MOUT,9002) CONV
      IF(EPH.GT.CONV) CONV=EPH
      WRITE(MOUT,9003) CONV
9020  J1=1
      SMALL=ONE/SMALL
      IF(IRANKB.EQ.0) GOTO 9036
      IF(SENS1.GE.ONE) GOTO 9035
      SENS1=SENS1*SMALL
      WRITE(MOUT,9004) J1,IRANKB,COND1,J1,IRANKB,SENS1
      GOTO 9036
9035  WRITE(MOUT,9049) J1,IRANKB,COND1,J1,IRANKB,SENS1,SMALL
9036  IF(IRANKB.EQ.IRANK) GOTO 9039
      J1=IRANKB+1
      IF(SENS2.GT.ONE) GOTO 9038
      SENS2=SENS2*SMALL
      WRITE(MOUT,9004) J1,IRANK,COND2,J1,IRANK,SENS2
      GOTO 9039
9038  WRITE(MOUT,9049) J1,IRANK,COND2,J1,IRANK,SENS2,SMALL
9039  WRITE(MOUT,9006) SIGDLH
      WRITE(MOUT,1001)
      IF(IFLAG.GT.0) WRITE(MOUT,9005)
      IF(IFLAG.LT.0) WRITE(MOUT,9905)
      DO 9031 K=1,M
9031  WRITE(MOUT,1003) T(K),(X(I,K),I=1,N)
      RETURN
C
C  FAIL EXIT MESSAGES
C---------------------
C
C  RANK-DEFICIENCY: BEST LEAST SQUARES SOLUTION OF BVP OBTAINED
9100  IFLAG=-1
      IF(KPRINT.LT.0) GOTO 9999
      WRITE(MOUT,9101)
      WRITE(MOUT,9002) CONVA
      IF(EPH.GT.CONVA) CONVA=EPH
      WRITE(MOUT,9003) CONVA
      IF(ITER.EQ.0) GOTO 9999
      SKAP=ZERO
      IF(FCA.EQ.ONE.AND.FC.EQ.ONE.AND.IRANKA.EQ.IRANK)
     @ SKAP=DSQRT(SUMXA/FCNUM)
      IF(SKAP.GT.ZERO) WRITE(MOUT,9102) SKAP
      GOTO 9999
C
C  TERMINATION AFTER MORE THAN ITMAX ITERATIONS
9200  IFLAG=-2
      IF(KPRINT.GE.0) WRITE(MOUT,9201) ITMAX
      GOTO 9999
C
C  SINGULAR TRAJECTORY
9310  WRITE(MOUT,9302)
9300  IFLAG=-3
      IF(KPRINT.LT.0) GOTO 9999
      J1=-KFLAG
      WRITE(MOUT,9301) J1,TFAIL
      GOTO 9999
C
C  CONVERGENCE FAIL OF GAUSS-NEWTON METHOD
9400  IFLAG=-4
      IF(KPRINT.GE.0) WRITE(MOUT,9401)
      GOTO 9999
C
C  INCONSISTENT INITIAL DATA
9500  IFLAG=-5
      IF(KPRINT.GE.0) WRITE(MOUT,9501)
      GOTO 9999
C
C  CONVERGENCE FAIL OF ITERATIVE REFINEMENT SWEEPS
9600  IFLAG=-6
      IF(KPRINT.LT.0) GOTO 9999
      WRITE(MOUT,9601)
      JN=JN-1
      IF(JN.GT.0) WRITE(MOUT,9602) JN
      GOTO 9999
C
C  INSUFFICIENT ERROR TOLERANCE FOR INTEGRATOR
9700  IFLAG=-7
      IF(KPRINT.LT.0) GOTO 9999
      TOLH=EPS/SIGDEL
      RELDIF=DSQRT(TOLH/SIGDEL)
      WRITE(MOUT,3203) TOLH,RELDIF
      WRITE(MOUT,9704) TOLH
      S=REDH/TOL
      DO 9710 J=1,M1
      IF(RF(J).GT.S) WRITE(MOUT,9705) J
9710  CONTINUE
      WRITE(MOUT,9703)
      GOTO 9999
C
C  ILL-CONDITIONED CONDENSED LINEAR SYSTEM
9800  IFLAG=-8
      IF(KPRINT.LT.0) GOTO 9999
      WRITE(MOUT,9801)
      IF(IRANK.EQ.NE) WRITE(MOUT,9802) EPX1H
      WRITE(MOUT,9803)
      GOTO 9999
C
C COMMON FAIL EXIT
9999  IF(KPRINT.LT.0) RETURN
      GOTO 9020
C
C
1001  FORMAT(1H1)
1002  FORMAT(14H0 INITIAL DATA,//)
1003  FORMAT(1H ,D13.5,3(D20.10),5(/14X,3(D20.10)))
1004  FORMAT(7H0 NODE=,I2,5H NBC=,I2,3H M=,I2,/,
     @   31H  PRESCRIBED RELATIVE PRECISION,
     @   D10.2,/,45H  MAXIMUM PERMITTED NUMBER OF ITERATION STEPS,I3//)
C
C  ITERATION MONITOR
1005  FORMAT(2H0 ,72(1H*))
1006  FORMAT(2H0 ,4X,2HIT,4X,2HNY,7X,6HLEVELF,9X,6HLEVELX,
     @                               7X,7HREL.FC.,3X,3HNEW,3X,6HRK-DEF)
2001  FORMAT(44H0 SINGULAR TRAJECTORY, RELAXATION FACTOR OR ,
     @                                      19HPSEUDO-RANK REDUCED,/)
3201  FORMAT(22H0 ITERATIVE REFINEMENT,/)
3202  FORMAT(7H0 SWEEP,I3,10H STARTS AT,I3,/,10(/,5D12.3))
3203  FORMAT(31H0 SUGGESTED INTEGRATOR ACCURACY,D10.1,/
     @      ,40H0 SUGGESTED RELATIVE DEVIATION PARAMETER,D10.1/)
3204  FORMAT(36H0 ADAPTED IN THE NEXT ITERATION STEP/)
3205  FORMAT(21H0 NORMS OF WRONSKIANS/10(5D12.3/))
4401  FORMAT(2H0 ,2(4X,I2),2(5X,D10.3),15X,I2,6X,I2)
5001  FORMAT(2H0 ,2(4X,I2),2(5X,D10.3),6X,F5.3)
6001  FORMAT(49H0 WARNING: ATTEMPT TO SOLVE INITIAL VALUE PROBLEM,/)
7701  FORMAT(2H0 ,4X,I2,31H NOT ACCEPTED RELAXATION FACTOR,
     @                                              11X,F5.3,12X,I2)
C
C  SOLUTION OUTPUT
9001  FORMAT(45H0 SOLUTION OF BOUNDARY VALUE PROBLEM OBTAINED,/,
     1 17H0 BVPLSQ REQUIRED,I3,21H ITERATION STEPS WITH,I4,
     2 23H TRAJECTORY EVALUATIONS,//)
9002  FORMAT(28H0 ACHIEVED RELATIVE ACCURACY,D10.3)
9003  FORMAT(28H0 RELIABLE RELATIVE ACCURACY,D10.3,/)
9004  FORMAT(16H0 SUBCONDITION (,I2,1H,,I2,2H) ,D10.3,/,
     @       16H0 SENSITIVITY  (,I2,1H,,I2,2H) ,D10.3,/)
9049  FORMAT(16H0 SUBCONDITION (,I2,1H,,I2,2H) ,D10.3,/,
     @       16H0 SENSITIVITY  (,I2,1H,,I2,2H) ,D10.3,2H *,1PD7.0/)
9005  FORMAT(15H0 SOLUTION DATA,/)
9006  FORMAT(29H0 MULTIPLE SHOOTING CONDITION,D10.3,/)
C
C  ERROR MESSAGES
9101  FORMAT(42H0 ITERATION TERMINATES AT STATIONARY POINT,/)
9102  FORMAT(30H0 INCOMPATIBILITY FACTOR KAPPA  ,D10.3,/)
9201  FORMAT(35H0 ITERATION TERMINATES AFTER ITMAX=,I3,
     @                                          17H  ITERATION STEPS)
9301  FORMAT(20H0 BVPLSQ TERMINATES ,/,
     @     13H0 SUBINTERVAL,I3,26H POSSIBLY INSERT NEW NODE ,D20.11/)
9302  FORMAT(52H0 SINGULAR TRAJECTORY BY DIFFERENCE APPROXIMATION OF,
     @                                     20H THE JACOBIAN MATRIX,/)
9401  FORMAT(39H0 GAUSS NEWTON METHOD FAILS TO CONVERGE,/)
9501  FORMAT(50H0 ERROR: INITIAL DATA AND BOUNDARY CONDITIONS ARE ,
     @                                             12HINCONSISTENT,/)
9601  FORMAT(46H0 TERMINATION SINCE ITERATIVE REFINEMENT FAILS,
     @                    12H TO CONVERGE,/,19H  INSERT NEW NODES ,/)
9602  FORMAT(17H  IN SUBINTERVAL ,I3,/)
9703  FORMAT(47H0 RELIABLE RELATIVE ACCURACY GREATER THAN 1.D-2,/)
9704  FORMAT(51H0 REDUCE RELATIVE ERROR TOLERANCE FOR INTEGRATOR TO,
     @            D10.1,/,21H  OR INSERT NEW NODES,/)
9705  FORMAT(17H  IN SUBINTERVAL ,I3,/)
9801  FORMAT(34H0 GAUSSIAN BLOCK ELIMINATION FAILS,/,
     @       44H  BY ILL-CONDITIONED CONDENSED LINEAR SYSTEM,/)
9802  FORMAT(26H0 RELATIVE ACCURACY OF DX1,D10.3)
9803  FORMAT(48H0 POSSIBLY TURN TO CODE BVPSOG INSTEAD OF BVPLSQ,/)
9905  FORMAT(12H0 FINAL DATA,/)
C
C  END SUBROUTINE BFBVPL
C
      END
C
      SUBROUTINE BFDERA (BC,N,NBC,M,XW,X1,XM,R,RH,A,B,RELDIF)
C-----------------------------------------------------------------------
C  DIFFERENCE APPROX. OF BOUNDARY DERIV. MATRICES A(NBC,N) AND B(NBC,N)
C-----------------------------------------------------------------------
      DOUBLE PRECISION  RELDIF,S,ONE,XH,ZERO
      DOUBLE PRECISION XW(N,M),A(NBC,N),B(NBC,N),R(NBC),RH(NBC),
     @                 X1(N),XM(N)
C
      EXTERNAL BC
C
      DATA  ONE/1.D0/ , ZERO/0.D0/
C
      DO 10 K=1,N
      XH=X1(K)
      S=RELDIF*XW(K,1)
      IF(XH.LT.ZERO) S=-S
      X1(K)=XH+S
      CALL BC (X1,XM,RH)
      X1(K)=XH
      S=ONE/S
      DO 102 I=1,NBC
102   A(I,K)=(RH(I)-R(I))*S
      XH=XM(K)
      S=RELDIF*XW(K,M)
      IF(XH.LT.ZERO) S=-S
      XM(K)=XH+S
      CALL BC (X1,XM,RH)
      XM(K)=XH
      S=ONE/S
      DO 103 I=1,NBC
103   B(I,K)=(RH(I)-R(I))*S
10    CONTINUE
C
      RETURN
C
C  END SUBROUTINE BFDERA
C
      END
      SUBROUTINE BFDERG (FCN,N,NE,M,M1,T,X,XU,XW,XJ,TJ,G,
     @                         ICOL,IVPSOL,HSTART,TOL,RELDIF,KFLAG)
C----------------------------------------------------------------------
C  DIFFERENCE APPROXIMATION OF WRONSKIAN MATRICES G(1),..,G(M1)
C----------------------------------------------------------------------
      DOUBLE PRECISION HSTART,HMAX,H,HSAVE,ONE,S,RELDIF,TJ,TJ1,TJA,TH
     @                ,TOL,DABS,ZERO
      DOUBLE PRECISION T(M),X(N,M), G(N,N,M1), XW(N,M), XU(N,M1), XJ(N)
      INTEGER ICOL(N)
C
      EXTERNAL FCN,IVPSOL
C
      DATA  ONE/1.D0/ , ZERO/0.D0/
C
      HSAVE=HSTART
      J=1
50    J1=J+1
      TJA=T(J)
      TJ1=T(J1)
      HMAX=DABS(TJ1-TJA)
      DO 500 IK=1,N
      I=ICOL(IK)
      H=HSAVE
      IF(J.EQ.1.AND.IK.GT.NE) GOTO 500
      TJ=TJA
      KFLAG=0
      DO 503 K=1,N
503   XJ(K)=X(K,J)
      TH=XJ(I)
      S=XW(I,J)*RELDIF
      IF(TH.LT.ZERO) S=-S
      XJ(I)=TH+S
      S=ONE/S
      CALL IVPSOL (N,FCN,TJ,XJ,TJ1,TOL,HMAX,H,KFLAG)
      IF(H.EQ.ZERO) GOTO 999
      DO 52 K=1,N
52    G(K,I,J)=S*(XJ(K)-XU(K,J))
500   CONTINUE
      HSAVE=H
      J=J1
      IF(J.LT.M) GOTO 50
C
      KFLAG=0
      RETURN
999   KFLAG=-J
C  ERROR RETURN
      RETURN
C
C  END SUBROUTINE BFDERG
C
      END
      SUBROUTINE BFRK1G (N,M,M1,XW,DX,HH,HHA,DXJ,G,FCA)
C-----------------------------------------------------------------------
C  RANK-1 UPDATES OF WRONSKIAN MATRICES G(1),...,G(M1)
C-----------------------------------------------------------------------
       DOUBLE PRECISION DNM,FCA,FCH,ONE,S,T,ZERO
       DOUBLE PRECISION G(N,N,M1), DX(N,M),XW(N,M), HH(N,M1),HHA(N,M1)
     1                 ,DXJ(N)
C
      DATA  ZERO/0.D0/ , ONE/1.D0/
C
      FCH=FCA-ONE
C
      DO 100 J=1,M1
      DNM=ZERO
      DO 110 I=1,N
      T=DX(I,J)/XW(I,J)
      DXJ(I)=T/XW(I,J)
110   DNM=DNM+T*T
      DNM=DNM*FCA
      IF(DNM.EQ.ZERO) GOTO 100
      DO 120 K=1,N
      T=DXJ(K)/DNM
      DO 120 I=1,N
      S=G(I,K,J)
      IF(S.NE.ZERO) G(I,K,J)=S+T*(HH(I,J)+FCH*HHA(I,J))
120   CONTINUE
100   CONTINUE
      RETURN
C
C  END SUBROUTINE BFRK1G
C
      END
      SUBROUTINE BFRHS1 (N,NBC,NRE,NE,M1,JIN,HH,R,B,G,U,DE,V,BG,IROW)
C----------------------------------------------------------------------
C  COMPUTATION OF CONDENSED RIGHT-HAND SIDE U(NRE)
C----------------------------------------------------------------------
      DOUBLE PRECISION S,TH,ZERO
      DOUBLE PRECISION G(N,N,M1), HH(N,M1), B(NBC,N), BG(NBC,N)
     @                ,U(NBC),DE(NBC),V(N),R(NBC)
      INTEGER IROW(N)
C
      DATA  ZERO/0.D0/
C
      DO 100 I=1,NRE
      IR=IROW(I)
100   U(I)=DE(IR)*R(IR)
      IF(JIN.GT.M1) RETURN
C
      DO 110 I=1,NRE
      IR=IROW(I)
      S=U(I)
      DO 111 K=1,N
      TH=DE(IR)*B(IR,K)
      BG(I,K)=TH
111   S=S+TH*HH(K,M1)
110   U(I)=S
      IF(M1.EQ.1.OR.JIN.EQ.M1) RETURN
C
      M2=M1-1
      DO 200 JJ=JIN,M2
      J=M2+JIN-JJ
      J1=J+1
      DO 200 I=1,NRE
      DO 210 K=1,N
      S=ZERO
      DO 211 L=1,N
211   S=S+BG(I,L)*G(L,K,J1)
210   V(K)=S
      S=U(I)
      DO 220 K=1,N
      S=S+V(K)*HH(K,J)
220   BG(I,K)=V(K)
200   U(I)=S
C
C  END SUBROUTINE BFRHS1
C
      RETURN
      END
      SUBROUTINE BFSCLE (N,M,M1,X,XU,XW,XTHR)
C----------------------------------------------------------------------
C  PROVIDES SCALING XW(N,M) OF VARIABLES X(N,M)
C----------------------------------------------------------------------
C
      DOUBLE PRECISION X(N,M),XW(N,M), XU(N,M1)
      DOUBLE PRECISION DABS,EPMACH,HALF,RED,SMALL,XMAX,XTHR,ZERO
C
      COMMON /MACHIN/ EPMACH, SMALL
C
      DATA ZERO/0.D0/, HALF/0.5D0/, RED/1.D-2/
C
      DO 220 I=1,N
220   XW(I,1)=DABS(X(I,1))
C
C  ARITHMETIC MEAN FOR XW(N,2),...,XW(N,M)
      DO 221 J=1,M1
      J1=J+1
      DO 221 I=1,N
221   XW(I,J1)=(DABS(X(I,J1))+DABS(XU(I,J)))*HALF
C
C  THRESHOLD
      DO 222 I=1,N
      XMAX=ZERO
      DO 223 J=1,M
      IF (XMAX.LT.XW(I,J)) XMAX=XW(I,J)
223   CONTINUE
      XMAX=XMAX*RED
      IF (XMAX.LT.XTHR) XMAX=XTHR
      DO 224 J=1,M
      IF(XW(I,J).LT.XMAX) XW(I,J)=XMAX
224   CONTINUE
222   CONTINUE
C
C  END SUBROUTINE BFSCLE
C
      RETURN
      END
      SUBROUTINE BFRECU (N,M,M1,JIN,HH,G,DX,U,V)
C---------------------------------------------------------------------
C  RECURSIVE SOLUTION OF M1 LINEAR (N,N)-SYSTEMS
C---------------------------------------------------------------------
      DOUBLE PRECISION S,ZERO
      DOUBLE PRECISION G(N,N,M1), DX(N,M), HH(N,M1), U(N),V(N)
C
C
      DATA  ZERO/0.D0/
C
      DO 10 I=1,N
10    U(I)=DX(I,1)
C
      DO 100 J=1,M1
      J1=J+1
      DO 110 I=1,N
      S=ZERO
      IF(J.GE.JIN) S=HH(I,J)
      DO 111 K=1,N
111   S=S+G(I,K,J)*U(K)
      V(I)=S
110   DX(I,J1)=S
      DO 120 I=1,N
120   U(I)=V(I)
100   CONTINUE
C
C  END SUBROUTINE BFRECU
C
      RETURN
      END
      SUBROUTINE BFPRJT (N,NE,IRANK,DEL1,U,D,V,QE,PIVOT)
C----------------------------------------------------------------------
C  TO BE USED IN CONNECTION WITH EITHER DECOMP/SOLVE OR BVDECC/BVSOLC
C----------------------------------------------------------------------
      DOUBLE PRECISION DEL1,S,ZERO, QE(N,N), U(N),D(N),V(N)
      INTEGER PIVOT(N)
C
      DATA  ZERO/0.D0/
C
      DO 10 L=1,NE
      I=PIVOT(L)
10    V(L)=U(I)
C
      IRK1=IRANK+1
      DEL1=ZERO
      DO 100 I=IRK1,NE
      S=V(I)
      DO 110 L=1,IRANK
110   S=S-V(L)*QE(L,I)
      IF(I.EQ.IRK1) GOTO 122
      I1=I-1
      DO 120 L=IRK1,I1
120   S=S-QE(L,I)*V(L)
122   S=S/D(I)
      DEL1=DEL1-S*S
100   V(I)=S
C
      DO 200 I=IRK1,NE
      K=NE+IRK1-I
      S=V(K)
      IF(K.EQ.NE) GOTO 212
      K1=K+1
      DO 210 L=K1,NE
210   S=S-QE(K,L)*V(L)
212   S=S/D(K)
200   V(K)=S
C
      DO 300 I=1,IRANK
      S=ZERO
      DO 310 L=IRK1,NE
310   S=S+V(L)*QE(I,L)
300   V(I)=-S
C
      DO 400 L=1,NE
      I=PIVOT(L)
400   U(I)=V(L)
C
      RETURN
C     END SUBROUTINE BFPRJT
C
      END
C
C*    Group  Initial value problem solver (Code BFDFX1)   
C
      SUBROUTINE BFDFX1 (N,F,X,Y,XEND,EPS,HMAX,H,KFLAG)
C
C
C  EXPLICIT EXTRAPOLATION INTEGRATOR
C  FOR NON-STIFF SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS
C  (BASED ON THE EXPLICIT MID-POINT DISCRETIZATION)
C
C**********  REVISION 1  *******  LATEST CHANGE: JUNE 4,82  ************
C
C  INTERNAL OPTION FOR POLYNOMIAL OR RATIONAL EXTRAPOLATION
C
C
C  REFERENCES:
C
C /1/ W.B.GRAGG:
C     ON EXTRAPOLATION ALGORITHMS FOR ORDINARY INITIAL VALUE PROBLEMS
C     SIAM J. NUMER. ANAL. 2, 384-404 (1965)
C
C /2/ R.BULIRSCH, J.STOER:
C     NUMERICAL TREATMENT OF ORDINARY DIFFERENTIAL EQUATIONS BY
C     EXTRAPOLATION METHODS
C     NUMER. MATH. 8, 1-13 (1966)
C
C /3/ P.DEUFLHARD:
C     ORDER AND STEPSIZE CONTROL IN EXTRAPOLATION METHODS
C     UNIVERSITY OF HEIDELBERG, SFB 123: TECH. REP. 93 (1980)
C
C
C  EXTERNAL SUBROUTINE (TO BE SUPPLIED BY THE USER)
C
C    F(X,Y,DY)          RIGHT-HAND SIDE OF FIRST-ORDER
C                       DIFFERENTIAL EQUATIONS
C      N  (.LE.13)      NUMBER OF FIRST-ORDER ODE'S
C      X                ACTUAL POSITION
C      Y(N)             VALUES AT X
C      DY(N)            DERIVATIVES AT X
C
C
C  INPUT PARAMETERS(* MARKS INOUT PARAMETERS)
C
C    N       (.LE.13)   NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS
C  * X                  STARTING POINT OF INTEGRATION
C  * Y(N)               INITIAL VALUES Y(1),...,Y(N)
C    XEND               PRESCRIBED FINAL POINT OF INTEGRATION
C    EPS                PRESCRIBED RELATIVE PRECISION (.GT.0)
C    HMAX               MAXIMUM PERMITTED STEPSIZE
C  * H                  INITIAL STEPSIZE GUESS
C  * KFLAG              PRINT PARAMETER
C                        0   NO OUTPUT
C                        1   INTEGRATION MONITOR
C                        2   ADDITIONALLY INTERMEDIATE SOLUTION POINTS
C                            T,Y(I),I=1,N
C
C                    (OUTPUT IS WRITTEN ON LOGICAL UNIT LOUT=6)
C
C  OUTPUT PARAMETERS
C
C    X                  ACHIEVED FINAL POINT OF INTEGRATION
C    Y(N)               FINAL VALUES
C    H                  STEPSIZE PROPOSAL FOR NEXT INTEGRATION STEP
C                       (H.EQ.0. ,IF DIFEX1 FAILS TO PROCEED)
C    KFLAG       .GE. 0:SUCCESSFUL INTEGRATION
C                       (KFLAG NOT ALTERED INTERNALLY)
C                .EQ.-1:MORE THAN NSTMAX BASIC INTEGRATION STEPS PER
C                       INTERVAL HAVE BEEN PERFORMED
C                .EQ.-2:MORE THAN JRMAX STEPSIZE REDUCTIONS
C                       OCCURRED PER BASIC INTEGRATION STEP
C                .EQ.-3:STEPSIZE PROPOSAL FOR NEXT BASIC INTEGRATION
C                       TOO SMALL
C
C
      INTEGER NJ(10),INCR(10),NRED( 9)
C     INTEGER NJ(JM),INCR(JM),NRED(KM)
      DOUBLE PRECISION Y(N),YL(13),YM(13),DY(13),DZ(13),S(13),DT(13,10),
C                      Y(N),YL(N) ,YM(N) ,DY(N) ,DZ( N),S( N),DT( N,JM),
     1  D(10,10),A(10),AL(10,10)
C       D(JM,JM),A(JM),AL(JM,JM)
      DOUBLE PRECISION B,B1,C,DABS,DBLE,DIFF,DM,DMAX,DSQRT,EPH,EPMACH,
     1 EPS,ERR,FC,FCM,FCO,FIVE,FJ,FJ1,FMIN,FN,G,H,HALF,HMAX,HMAXU,HR,
     2 H1,OMJ,OMJO,ONE,ONE1,Q,QUART,RED,RO,SAFE,SK,SMALL,TA,TEN,
     3 U,V,W,X,XEND,XEPS,XN,YH,ZERO
C
      COMMON /COUNT/NSTEP,NFCN
C
      DATA LOUT/6/
C
      DATA ZERO/0.D0/,FMIN/1.D-2/,RO/0.25D0/,QUART/0.25D0/,HALF/0.5D0/,
     *     SAFE/0.5D0/,ONE/1.D0/,ONE1/1.01D0/,TEN/1.D1/,
     *     FIVE/5.D0/
      DATA DT/130*0.D0/
C     DATA DT/N*JM*0.D0/
C
C  STEPSIZE SEQUENCE HA (DUE TO /3/ )
      DATA NJ/2,4,6,8,10,12,14,16,18,20/
C
      EXTERNAL F
C
C  RELATIVE MACHINE PRECISION
      CALL ZIBCONST(EPMACH,SMALL)
C
C  ASSOCIATED MAXIMUM COLUMN NUMBER (1.LE.KM.LE.9)
      KM=6
C
C  ASSOCIATED MAXIMUM ROW NUMBER (2.LE.JM.LE.10)
      JM=KM+1
C
C  INTERNAL PARAMETERS
C  STANDARD VALUES FIXED BELOW
C
C
C
C  POLYNOMIAL EXTRAPOLATION (AITKEN-NEVILLE ALGORITHM)
C  (RECOMMENDED STANDARD OPTION)
      IPOL=1
C  RATIONAL EXTRAPOLATION (STOER ALGORITHM)
C     IPOL=0
C
C  MAXIMUM PERMITTED NUMBER OF INTEGRATION STEPS PER INTERVAL
      NSTMAX=50000
C
C  MAXIMUM PERMITTED NUMBER OF STEPSIZE REDUCTIONS
      JRMAX=5
C
C  INITIAL PREPARATIONS
      EPH=RO*EPS
      FJ1=DBLE(NJ(1))
      A(1)=FJ1+ONE
      DO 60 J=2,JM
      J1=J-1
      INCR(J1)=0
      NRED(J1)=0
      FJ=DBLE(NJ(J))
      V=A(J1)+FJ
      A(J)=V
      DO 61 K=1,J1
      W=FJ/DBLE(NJ(K))
61    D(J,K)=W*W
      IF(J.EQ.2) GOTO 60
      W=V-FJ1
      DO 62 K1=2,J1
      K=K1-1
      U=(A(K1)-V)/(W*DBLE(K+K1))
      U=EPH**U
62    AL(J1,K)=U
60    CONTINUE
      KOH=1
      JOH=2
65    IF(JOH.GE.JM) GOTO 66
      IF(A(JOH+1)*ONE1.GT.A(JOH)*AL(JOH,KOH)) GOTO 66
      KOH=JOH
      JOH=JOH+1
      GOTO 65
66    K=0
      KM=KOH
      JM=KM+1
      INCR(JM)=-1
      OMJO=ZERO
      IF(KFLAG.GE.1) WRITE(LOUT,1001) EPS,KM
      EPMACH=EPMACH*TEN
      HMAX=DABS(HMAX)
      NSTEP=0
      NFCN=0
      XEPS=(DABS(X)+DABS(XEND))*EPMACH
      FN=DBLE(N)
      H1=XEND-X
      HMAXU=HMAX
      HR=HMAX
      DMAX=FIVE
C
C  BASIC INTEGRATION STEP
401   IF(DABS(H1).LE.XEPS) GOTO 403
      Q=H1/H
      IF(Q.LE.EPMACH) GOTO 403
      IF(KFLAG.GT.1) WRITE(LOUT,1009) NSTEP,NFCN,X,K,KOH
      IF(KFLAG.GT.1) WRITE(LOUT,1000) NSTEP,NFCN,X,(Y(I),I=1,N)
      IF(Q.GE.ONE1) GOTO 402
      HR=H
      H=H1
402   JRED=0
      DO 405 K=1,KM
405   INCR(K)=INCR(K)+1
      HMAX=DABS(H1)
      IF(HMAXU.LT.HMAX) HMAX=HMAXU
C
C  SCALING
C  (FOR REAL LIFE APPLICATIONS TO BE POSSIBLY ALTERED BY THE USER)
      DO 5 I=1,N
      U=DABS(Y(I))
      IF(U.LT.EPMACH) U=ONE
5     S(I)=U
C
C  EXPLICIT EULER STARTING STEP
      CALL  F (X,Y,DZ)
      NFCN=NFCN+1
C
10    XN=X+H
      FCM=DABS(H)/HMAX
      IF(FCM.LT.FMIN) FCM=FMIN
C
      DO 260 J=1,JM
      M=NJ(J)
      G=H/DBLE(M)
      B=G+G
      DO 210 I=1,N
      YL(I)=Y(I)
210   YM(I)=Y(I)+G*DZ(I)
      M=M-1
C  EXPLICIT MID-POINT RULE
      DO 220 K=1,M
      CALL  F (X+G*DBLE(K),YM,DY)
      NFCN=NFCN+1
      DO 220 I=1,N
      U=YL(I)+B*DY(I)
      YL(I)=YM(I)
      YM(I)=U
220   CONTINUE
C  FINAL STEP
      CALL  F (XN,YM,DY)
      NFCN=NFCN+1
      DM=ZERO
      DO 2200 I=1,N
      YH=YL(I)+G*DY(I)
      DIFF=YH-YM(I)
      YM(I)=(YM(I)+YH)*HALF
      DIFF=DABS(DIFF)
      SK=DABS(YM(I))
      IF(SK.LT.S(I)) SK=S(I)
      DIFF=DIFF/SK
      IF(DIFF.GT.DM) DM=DIFF
2200  CONTINUE
C
C  STABILITY CHECK
      IF(DM.LT.DMAX) GOTO 2209
C
C  EMERGENCY EXIT
      IF(KFLAG.GT.0) WRITE(LOUT,1006)
      GOTO 2601
C
C  PREVENTION OF POSSIBLE ORDER INCREASE
2209  IF(J.GT.2.OR.DM.LT.DMAX*HALF) GOTO 2207
      DO 2208 L=JOH,JM
      IF(INCR(L).GT.0) INCR(L)=0
      INCR(L)=INCR(L)-2
2208  CONTINUE
C
C  EXTRAPOLATION
C  ( IPOL=1: POLYNOMIAL, IPOL=0: RATIONAL)
2207  ERR=ZERO
      DO 234 I=1,N
      V=DT(I,1)
      C=YM(I)
      DT(I,1)=C
      IF(J.EQ.1) GOTO 234
      TA=C
      DO 231 K=2,J
      JK=J-K+1
      B1=D(J,JK)
      W=C-V
      IF(IPOL.EQ.0) GOTO 229
      U=W/(B1-ONE)
      C=B1*U
      GOTO 230
229   B1=B1*V
      B=B1-C
      U=V
      IF(B.EQ.ZERO) GOTO 230
      B=W/B
      U=C*B
      C=B1*B
230   V=DT(I,K)
      DT(I,K)=U
231   TA=U+TA
      YM(I)=TA
      TA=DABS(TA)
      IF(TA.LT.S(I)) TA=S(I)
      U=U/TA
      ERR=ERR+U*U
234   CONTINUE
      IF(J.EQ.1) GOTO 260
C ERROR (SCALED ROOT MEAN SQUARE)
      ERR=DSQRT(ERR/FN)
      KONV=0
      IF(ERR.LT.EPS) KONV=1
      ERR=ERR/EPH
C
C ORDER CONTROL
      K=J-1
      L=J+K
      FC=ERR**(ONE/DBLE(L))
      IF(FC.LT.FCM) FC=FCM
C  OPTIMAL ORDER DETERMINATION
      OMJ=FC*A(J)
      IF(J.GT.2.AND.OMJ*ONE1.GT.OMJO.OR.K.GT.JOH) GOTO 235
      KO=K
      JO=J
      OMJO=OMJ
      FCO=FC
235   IF(J.LT.KOH.AND.NSTEP.GT.0) GOTO 260
      IF(KONV.EQ.0) GOTO 236
      IF(KO.LT.K.OR.INCR(J).LT.0) GOTO 20
C  POSSIBLE INCREASE OF ORDER
      IF(NRED(KO).GT.0) NRED(KO)=NRED(KO)-1
      FC=FCO/AL(J,K)
      IF(FC.LT.FCM) FC=FCM
      J1=J+1
      IF(A(J1)*FC*ONE1.GT.OMJO) GOTO 20
      FCO=FC
      KO=JO
      JO=JO+1
      GOTO 20
C
C
C  CONVERGENCE MONITOR
236   RED=ONE/FCO
      JK=KM
      IF(JOH.LT.KM) JK=JOH
      IF(K.GE.JK) GOTO 239
      IF(KO.LT.KOH) RED=AL(KOH,KO)/FCO
237   IF(AL(JK,KO).LT.FCO) GOTO 239
260   CONTINUE
C
C STEPSIZE REDUCTION (DUE TO EXTRAPOLATION TABLE)
239   RED=RED*SAFE
      H=H*RED
2392  IF(NSTEP.EQ.0) GOTO 2390
      NRED(KOH)=NRED(KOH)+1
      DO 2391 L=KOH,KM
2391  INCR(L)=-2-NRED(KOH)
2390  JRED=JRED+1
      IF(KFLAG.GT.0) WRITE(LOUT,1002) JRED,RED,KOH
      IF(JRED.GT.JRMAX) GOTO 32
      GOTO 10
C
C  STEPSIZE REDUCTION (DUE TO STABILITY)
2601  HMAX=G*FJ1*QUART
      RED=HMAX/DABS(H)
      H=HMAX
      IF(JRED.GT.0) GOTO 2390
      GOTO 2392
C
C  PREPARATIONS FOR NEXT BASIC INTEGRATION STEP
20    X=XN
      H1=XEND-X
      DO  2606 I=1,N
2606  Y(I)=YM(I)
      NSTEP=NSTEP+1
      IF(NSTEP.GT.NSTMAX) GO TO 31
C
C STEPSIZE PREDICTION
      IF(FCO.NE.FCM) HR=H
      H=H/FCO
      KOH=KO
      JOH=KOH+1
      IF(DABS(H).GT.DABS(X)*EPMACH) GO TO 401
      GO TO 33
C
C  SOLUTION EXIT
403   H=HR
      IF(KFLAG.GT.1) WRITE(LOUT,1009) NSTEP,NFCN,X,K,KOH
      IF(KFLAG.GT.1) WRITE(LOUT,1000) NSTEP,NFCN,X,(Y(I),I=1,N)
      HMAX=HMAXU
      RETURN
C
C  FAIL EXIT
31    IF(KFLAG.GE.1)WRITE(LOUT,1008) NSTMAX
      KFLAG=-1
      GOTO 39
32    IF(KFLAG.GE.1) WRITE(LOUT,1010)JRMAX
      KFLAG=-2
      GOTO 39
33    IF(KFLAG.GE.1) WRITE(LOUT,1004)
      KFLAG=-3
39    H=ZERO
      HMAX=HMAXU
      RETURN
C
1000  FORMAT(1H ,2I9,5D20.11,/,(1H ,38X,4D20.11))
1001  FORMAT(1H0,25HDIFEX1     REL.PREC. EPS ,D10.3,8HMAX.COL.,I3
     *        ,//,5X,4HSTEP,3X,7HF-CALLS,8X,1HX,19X,7HY1(X)..,//)
1002  FORMAT(1H ,I3,17HREDUCTION FACTOR ,D10.3,I9,/)
1004  FORMAT(/,40H0  STEPSIZE REDUCTION FAILED TO SUCCEED, //)
1006  FORMAT(/,46H0  STABILITY CHECK ACTIVATED  STEPSIZE REDUCED,/)
1008  FORMAT(18H0MORE THAN NSTMAX=,I6,18H INTEGRATION STEPS,//)
1009  FORMAT(1H ,2I9,D20.11,I9,I6,/)
1010  FORMAT(17H0MORE THAN JRMAX=,I3,29H STEPSIZE REDUCTIONS PER STEP,/)
C
C
C
C  END BFDFX1
C
      END
