PROGRAM ABCALL
      PARAMETER(MAXDEG=13,MXDM=51)
C N, THE NUMBER OF EQUATIONS IN THE SYSTEM F(U,R) = 0
C MAXDEG, MAXIMUM DEGREE PREDICTOR FOR WHICH STORAGE IS ALLOCATED
C MXDM, THE LEADING DIMENSION OF TWO DIMENSIONAL ARRAYS D,PHI,PHIA
C       STORAGE IS ALLOCATED FOR SYSTEMS OF UP TO MXDM-1 EQUATIONS
C********************************************************************
C
C ** TYPE AND DIMENSION**
      REAL WC(MXDM),WCP(MXDM),WP(MXDM),DIR(MXDM),U(MXDM),FFIN(MXDM),
     +    V(MXDM),STO(MXDM),STO2(MXDM),D(MXDM,MXDM),TANGNT(MXDM), 
C ARRAYS FOR PREDICTOR AND SCALED DIVIDED DIFFERENCE STORAGE
     +     G(0:MAXDEG,1:MAXDEG),PSYP(0:MAXDEG),PSY(0:MAXDEG),
     +     A(0:MAXDEG),B(0:MAXDEG),PHI(MXDM,0:MAXDEG),
     +     PSYARC(0:MAXDEG),PHIA(MXDM,0:MAXDEG),
C REAL SPECIFICATION VARIABLES
     +     EREL,EABS,ECOR,EFIN,HMIN,HMAX,RLOW,RHIGH,ORIENT,DS
C
      INTEGER IPIVOT(MXDM),IFLAG,KNTROL,LEX,MXDG,MXSTEP,
     +        NPRNT,M,NEWTON,NFAC,IFIN,NP1,MPREF,COUNT,ITOK
C
      EXTERNAL F1,DZF1,DRF1
      EXTERNAL F2,DZF2,DRF2
      EXTERNAL F3,DZF3,DRF3
      EXTERNAL F4,DZF4,DRF4
      EXTERNAL F5,DZF5,DRF5
      EXTERNAL F6,DZF6,DRF6
      EXTERNAL F7,DZF7,DRF7
      EXTERNAL F8,DZF8,DRF8
C
C  INITIALIZE PARAMETERS AND INPUT VARIABLES
C *********** PROBLEM 1 WATSTON N=10 ***************************** 
      PRINT*,'PROBLEM 1 WATSON N = 10'
      N=10
C ENTER INITIAL SOLUTION POINT TO F(W) = 0.
      DO 29 L=1,N+1
         WC(L) = 0.E0
29    CONTINUE
C** PARAMETERS ********************
C THE BOUNDS OF CONTINUATION PARAMETER: THE CONTINUATION WILL
C STOP WITH A POINT WC WITH WC(N+1) = RHIGH, (OR RLOW).
      RHIGH = 1.
      RLOW = -20.
C MINIMUM PREDICTOR STEPSIZE
      HMIN = 1.E-6
C MAXIMUM PREDICTOR STEPSIZE
      HMAX = 100.
C ABSOLUTE AND RELATIVE PREDICTOR ERROR PARAMETERS
      EABS = 1.E-2
      EREL = 1.E-2
C ERROR PARAMETER FOR CORRECTED POINTS
      ECOR = 1.E-4
C ERROR PARAMETER FOR FINAL TARGET POINT:
C IF EFIN .GT. 0, USE A CHORD CORRECTOR, STOPPING ERROR EFIN
C IF EFIN .LT. 0, USE A NEWTON CORRECTOR, STOPPING ERROR -EFIN
C IF EFIN .EQ. 0, USE NO CORRECTOR
      EFIN = +1.0E-4
C CONTROL PARAMETER FOR FINAL POINT CORRECTION
C IFIN .GT. 0,  FINAL CORRECTION DONE WITH A NEWTON OR CHORD ITERATION
C             WITH TARGET PARAMETER FIXED AND NO BORDERING AND .LE IFIN
C               ITERATIONS.  
C IFIN .LT. 0, USUAL BORDERING CORRECTOR USED BUT WITH DIR RESET
C              TO GIVE CORRECTION WITH TARGET PARAMETER FIXED
C              AND .LE. -IFIN ITERATIONS
C IFIN .EQ 0,  USUAL CORRECTION FOR POINTS; BORDERING ALGORITHM
C              WITH ITERATES IN A HYP. PLANE ORTH. TO USUAL PREDICTOR.
C              MAXIMUM ITERATES DETERMINED AS WITH POINTS (SEE NEWTON)
      IFIN = -1
C SET MAXIMUM ITERATIONS IN CORRECTOR:
C IF NEWTON = 1, 2, ...  CORRECTOR ABORTS AFTER NEWTON CHORD ITERATES
C IF NEWTON = -1, -2, ... CORRECTOR USES UP TO "NEWTON" NEWTON ITERATES
      NEWTON = 7
C SET MAX CORRECTOR ITERATES FOR INCREASING STEPSIZE AFTER FAILURES
      ITOK = 7
C SET PRINT CODE:
C IF NPRNT=0, NO PRINTING OF INFORMATION
C IF NPRNT=1, INITIAL HEADING WITH PARAMETER SETTINGS PRINTED
C             FIRST TWO STEPS PRINTED ALONG WITH SINGULAR POINT
C             INFO AND FINAL POINT INFO
C IF NPRNT=2, ABOVE PLUS CORRECTED POINTS (FIRST TWO, EVERY KSKIP STEP)
C IF NPRNT=3, ABOVE PLUS HEADING FOR STEPS AND OTHER STEP INFO.
      NPRNT = -3
C SET SKIPPING PARAMETER: PRINT EVERY KSKIP STEP.
      KSKIP = 0
C SET THE MINIMUM NUMBER OF STEPS BEFORE FINAL POINT STOPPING
C WILL BE INVOKED.  THIS CAN BE USED TO PICK OUT SUBSEQUENT
C PASSINGS OF A PARAMETER VALUE FOR ATTENTION.  
C IF KTARG IS A NEGATIVE INTEGER, FINAL PREDICTED POINT RETURNED IN WC
      KTARG = 1
C SET ERROR TOLERANCE FOR MAX NORM OF F(W) AT INITIAL POINT
C IF THIS TOLERANCE IS NOT MET, THE CODE WILL ABORT.
      TOLINT = EABS
C SET ORIENTATION FOR INITIAL DIRECTION OF CONTINUATION
C IF ORIENT = +1, THE PARAMETER WC(N+1) INITIALLY IS INCREASING
C IF ORIENT = -1, THE PARAMETER INITIALLY DECREASES
      ORIENT = +1.
C SET THE INDEX OF THE TARGET COMPONENT FOR STOPPING THE CONTINUATION:
C IF NCOMP = 1, ..., N+1, THE CONTINUATION WILL STOP AND LOCATE A POINT
C            WITH WC(NCOMP) = VCOMP ONCE WC(NCOMP) .GT. VCOMP
C IF NCOMP = -1, ...,-(N+1) THE CONTINUATION WILL STOP AND LOCATE A
C            POINT WITH WC(-NCOMP) = VCOMP ONCE WC(-NCOMP) .LT. VCOMP
C IF NCOMP = 0, THIS FEATURE IS DISABLED (SEE RLOW, RHIGH HOWEVER)
      NCOMP = N+1
C SET TARGET VALUE FOR WC(ABS(MCOMP)) (SEE NCOMP)
      VCOMP = 1.0
C SET MAXIMUM NUMBER OF CONTINUATION STEPS (POINTS)
      MXSTEP = 600
C AUTOMATIC PREDICTOR ORDER INCREASE UNTIL ORDER > MPREF
      MPREF = 0
C  MAXIMUM DEGREE OF PREDICTOR POLYNOMIAL
      MXDG = 3
C NFAC = 0 FOR TANGENT COMPUTED AT PREDICTED POINT
C NFAC = 1 FOR TANGENT COMPUTED AT CORRECTED POINT
      NFAC = 0
C INITIAL STEPSIZE.  IF DS .LT.0, DS IS COMPUTED BY ABCON
      DS =  -8.E-2
C ISING = 0 TURNS OFF BIFURCATION CHECKING (FOLDS ONLY)
C ISING = 1 TURNS ON BUFURCATION POINT CHECKING
      ISING = 0
C
      CALL ABCON(N,WC,FFIN,F8,DZF8,DRF8,NPRNT,KSKIP,ORIENT,NCOMP,
     + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX,
     + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP,
     + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE)
C
CP      print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1)
      print*,'====================================================='
C######################################################################
C      GO TO 1000
C *********** PROBLEM 2 WATSON CURVE N = 12 ************** 
      PRINT*,'*********** PROBLEM 2 WATSON CURVE N = 12 **********'
      N=12
C ENTER INITIAL SOLUTION POINT TO F(W) = 0.
      DO 229 L=1,N+1
         WC(L) = 0.E0
229    CONTINUE
C** PARAMETERS ********************
      RHIGH = 1.
      RLOW = -20.
      HMIN = 1.E-6
      HMAX = 100.
      EABS = 1.E-2
      EREL = 1.E-2
      ECOR = 1.E-4
      EFIN = +1.0E-4
      IFIN = -1
      NEWTON = 7
      ITOK = 7
      NPRNT = -3
      KSKIP = 0
      KTARG = 1
      TOLINT = EABS
      ORIENT = +1.
      NCOMP = N+1
      VCOMP = 1.0
      MXSTEP = 800
      MPREF = 1
C     MXDG = 3
      NFAC = 0
      DS =  -8.E-2
      ISING = 0
C
      CALL ABCON(N,WC,FFIN,F8,DZF8,DRF8,NPRNT,KSKIP,ORIENT,NCOMP,
     + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX,
     + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP,
     + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE)
C
CP      print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1)
        print*,'====================================================='
C######################################################################
C *********** PROBLEM 3 (SCHWETLICK 3) N = 4 ************** 
      PRINT*,'****** PROBLEM 3 (SCHWETLICK 3) N = 4 ********'
      N=4
C ENTER INITIAL SOLUTION POINT TO F(W) = 0.
      DO 329 L=1,N+1
         WC(L) = 0.E0
329    CONTINUE
      WC(1) = -3.E0
      WC(2) = -1.E0
      WC(3) = -3.E0
      WC(4) = -1.E0
C** PARAMETERS ********************
      RHIGH = 1.
      RLOW = -20.
      HMIN = 1.E-8
      HMAX = 100.
      EABS = 1.E-5
      EREL = 1.E-5
      ECOR = 1.E-7
      EFIN = +1.0E-7
      IFIN = -1
      NEWTON = 7
      ITOK = 7
      NPRNT = -3
      KSKIP = 0
      KTARG = 1
      TOLINT = EABS
      ORIENT = +1.
      NCOMP = N+1
      VCOMP = 1.0
      MXSTEP = 300
      MPREF = 0
C     MXDG = 3
      NFAC = 0
      DS =  -8.E-2
C
      CALL ABCON(N,WC,FFIN,F1,DZF1,DRF1,NPRNT,KSKIP,ORIENT,NCOMP,
     + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX,
     + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP,
     + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE)
C
CP      print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1)
      print*,'====================================================='
C######################################################################
C *********** PROBLEM 4 (SCHWETLICK 4) N = 6 ************** 
      PRINT*,'********* PROBLEM 4 (SCHWETLICK 4) N = 6 ********'
      N=6
C ENTER INITIAL SOLUTION POINT TO F(W) = 0.
      DO 429 L=1,N+1
         WC(L) = 0.E0
429    CONTINUE
C** PARAMETERS ********************
      RHIGH = 1.
      RLOW = -20.
      HMIN = 1.E-8
      HMAX = 100.
      EABS = 1.E-3
      EREL = 1.E-3
      ECOR = 1.E-5
      EFIN = +1.0E-5
      IFIN = -1
      NEWTON = 7
      ITOK = 7
      NPRNT = -3
      KSKIP = 0
      KTARG = 52
      TOLINT = EABS
      ORIENT = +1.
      NCOMP = N+1
      VCOMP = 1.0
      MXSTEP = 100
      MPREF = 0
C     MXDG = 3
      NFAC = 0
      DS =  -8.E-2
C
      CALL ABCON(N,WC,FFIN,F2,DZF2,DRF2,NPRNT,KSKIP,ORIENT,NCOMP,
     + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX,
     + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP,
     + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE)
C
CP      print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1)
      print*,'====================================================='
C######################################################################
C *********** PROBLEM 5 (WATSON 1) N = 10 ************** 
      PRINT*,'********** PROBLEM 5 (WATSON 1) N = 10 ************'
      N=10
C ENTER INITIAL SOLUTION POINT TO F(W) = 0.
      DO 529 L=1,N+1
         WC(L) = 0.E0
529    CONTINUE
C** PARAMETERS ********************
      RHIGH = 1.
      RLOW = -20.
      HMIN = 1.E-6
      HMAX = 100.
      EABS = 1.E-2
      EREL = 1.E-2
      ECOR = 1.E-4
      EFIN = +1.0E-4
      IFIN = -1
      NEWTON = 7
      ITOK = 7
      NPRNT = -3
      KSKIP = 0
      KTARG = 1
      TOLINT = EABS
      ORIENT = +1.
      NCOMP = N+1
      VCOMP = 1.0
      MXSTEP = 20
      MPREF = 1
C     MXDG = 3
      NFAC = 0
      DS =  -8.E-2
C
      CALL ABCON(N,WC,FFIN,F4,DZF4,DRF4,NPRNT,KSKIP,ORIENT,NCOMP,
     + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX,
     + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP,
     + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE)
C
CP      print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1)
      print*,'====================================================='
C######################################################################
C *********** PROBLEM 6  N = 10 ************** 
      PRINT*,'*********** PROBLEM 6  N = 10 **************'
      N=10
C ENTER INITIAL SOLUTION POINT TO F(W) = 0.
      DO 629 L=1,N+1
         WC(L) = 0.E0
629    CONTINUE
C** PARAMETERS ********************
      RHIGH = 1.
      RLOW = -20.
      HMIN = 1.E-6
      HMAX = 100.
      EABS = 1.E-2
      EREL = 1.E-2
      ECOR = 1.E-4
      EFIN = +1.0E-4
      IFIN = -1
      NEWTON = 7
      ITOK = 7
      NPRNT = -3
      KSKIP = 0
      KTARG = 1
      TOLINT = EABS
      ORIENT = +1.
      NCOMP = N+1
      VCOMP = 1.0
      MXSTEP = 20
      MPREF = 1
C     MXDG = 3
      NFAC = 0
      DS =  -8.E-2
C
      CALL ABCON(N,WC,FFIN,F5,DZF5,DRF5,NPRNT,KSKIP,ORIENT,NCOMP,
     + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX,
     + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP,
     + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE)
C
CP      print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1)
      print*,'====================================================='
C######################################################################
C *********** PROBLEM 7 (BROWN FUNCTION) N = 10 ************** 
      PRINT*,'********* PROBLEM 7 (BROWN FUNCTION) N = 10 ********'
      N=10
C ENTER INITIAL SOLUTION POINT TO F(W) = 0.
      DO 729 L=1,N+1
         WC(L) = 0.E0
729    CONTINUE
C** PARAMETERS ********************
      RHIGH = 1.
      RLOW = -20.
      HMIN = 1.E-6
      HMAX = 100.
      EABS = 1.E-2
      EREL = 1.E-2
      ECOR = 1.E-4
      EFIN = +1.0E-4
      IFIN = -1
      NEWTON = 7
      ITOK = 7
      NPRNT = -3
      KSKIP = 0
      KTARG = 1
      TOLINT = EABS
      ORIENT = +1.
      NCOMP = N+1
      VCOMP = 1.0
      MXSTEP = 20
      MPREF = 1
C     MXDG = 3
      NFAC = 0
      DS =  -8.E-2
C
      CALL ABCON(N,WC,FFIN,F6,DZF6,DRF6,NPRNT,KSKIP,ORIENT,NCOMP,
     + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX,
     + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP,
     + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE)
C
CP      print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1)
      print*,'====================================================='
C######################################################################
C *********** PROBLEM 8 (BROWN FUNCTION) N = 25 ************** 
      PRINT*,'********* PROBLEM 8 (BROWN FUNCTION) N = 25 ********'
      N=25
C ENTER INITIAL SOLUTION POINT TO F(W) = 0.
      DO 829 L=1,N+1
         WC(L) = 0.E0
829    CONTINUE
C** PARAMETERS ********************
      RHIGH = 1.
      RLOW = -20.
      HMIN = 1.E-6
      HMAX = 100.
      EABS = 1.E-2
      EREL = 1.E-2
      ECOR = 1.E-4
      EFIN = +1.0E-4
      IFIN = -1
      NEWTON = 7
      ITOK = 7
      NPRNT = -3
      KSKIP = 0
      KTARG = 1
      TOLINT = EABS
      ORIENT = +1.
      NCOMP = N+1
      VCOMP = 1.0
      MXSTEP = 20
      MPREF = 1
C     MXDG = 3
      NFAC = 0
      DS =  -8.E-2
C
      CALL ABCON(N,WC,FFIN,F6,DZF6,DRF6,NPRNT,KSKIP,ORIENT,NCOMP,
     + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX,
     + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP,
     + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE)
C
CP      print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1)
      print*,'====================================================='
C ################################################################
C *********** PROBLEM 9 (BROWN FUNCTION) N = 50 ************** 
      PRINT*,'********* PROBLEM 9 (BROWN FUNCTION) N = 50 ********'
      N=50
C ENTER INITIAL SOLUTION POINT TO F(W) = 0.
      DO 929 L=1,N+1
         WC(L) = 0.E0
929    CONTINUE
C** PARAMETERS ********************
      RHIGH = 1.
      RLOW = -20.
      HMIN = 1.E-6
      HMAX = 100.
      EABS = 1.E-2
      EREL = 1.E-2
      ECOR = 1.E-4
      EFIN = +1.0E-4
      IFIN = -1
      NEWTON = 7
      ITOK = 7
      NPRNT = -3
      KSKIP = 0
      KTARG = 1
      TOLINT = EABS
      ORIENT = +1.
      NCOMP = N+1
      VCOMP = 1.0
      MXSTEP = 20
      MPREF = 1
C     MXDG = 3
      NFAC = 0
      DS =  -8.E-2
C
      CALL ABCON(N,WC,FFIN,F6,DZF6,DRF6,NPRNT,KSKIP,ORIENT,NCOMP,
     + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX,
     + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP,
     + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE)
C
CP      print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1)
      print*,'====================================================='
C###################################################################
C *********** PROBLEM 10 F AND R, REGULARIZING HOMOTOPY ************** 
      PRINT*,'PROBLEM 10 F AND R, REGULARIZING HOMOTOPY'
      N=2
C ENTER INITIAL SOLUTION POINT TO F(W) = 0.
      DO 1029 L=1,N+1
         WC(L) = 0.E0
1029    CONTINUE
      WC(1) = 15.
      WC(2) = -2.
      WC(3) = 0.
C** PARAMETERS ********************
      RHIGH = 1.
      RLOW = -20.
      HMIN = 1.E-6
      HMAX = 100.
      EABS = 1.E-2
      EREL = 1.E-2
      ECOR = 1.E-4
      EFIN = +1.0E-4
      IFIN = -1
      NEWTON = 7
      ITOK = 7
      NPRNT = -3
      KSKIP = 0
      KTARG = 1
      TOLINT = EABS
      ORIENT = +1.
      NCOMP = N+1
      VCOMP = 1.0
      MXSTEP = 100
      MPREF = 1
C     MXDG = 3
      NFAC = 0
      DS =  -8.E-2
C
      CALL ABCON(N,WC,FFIN,F7,DZF7,DRF7,NPRNT,KSKIP,ORIENT,NCOMP,
     + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX,
     + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP,
     + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE)
C
CP      print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1)
      print*,'====================================================='
C######################################################################
C *********** PROBLEM 11 F AND R, STANDARD HOMOTOPY ************** 
      PRINT*,'PROBLEM 11 F AND R, STANDARD HOMOTOPY'
C     do 999 l=1,5
C     NPRNT = l-1
C     print*,'nprnt =',nprnt
      N=2
C ENTER INITIAL SOLUTION POINT TO F(W) = 0.
      WC(1) = 15.
      WC(2) = -2.
      WC(3) = 0.
C  PARAMETERS
      RHIGH = 1.
      RLOW = -20.
      HMIN = 1.E-6
      HMAX = 100.
      EABS = 1.E-2
      EREL = 1.E-2
      ECOR = 1.E-4
      EFIN = +1.0E-4
      IFIN = -1
      NEWTON = 7
      ITOK = 7
      NPRNT = -3
      KSKIP = 0
      KTARG = 1
      TOLINT = EABS
      ORIENT = +1.
      NCOMP = N+1
      VCOMP = 1.0
      MXSTEP = 100
      MPREF = 0
      MXDG = 3
      NFAC = 0
      DS =  -8.E-2
C
      CALL ABCON(N,WC,FFIN,F3,DZF3,DRF3,NPRNT,KSKIP,ORIENT,NCOMP,
     + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX,
     + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP,
     + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE)
C
CP      print*,'IRCODE =',IRCODE,'wc=',(wc(i),i=1,n+1)
      print*,'====================================================='
999   continue
C OUTPUT RESULTS
1000   CONTINUE
C
      END
C      ABCON    BY    B. N. LUNDBERG AND A. B. POORE
C
      SUBROUTINE ABCON(N,WC,FFIN,F,DZF,DRF,NPRNT,KSKIP,ORIENT,NCOMP,
     + VCOMP,RHIGH,RLOW,KTARG,IFIN,EABS,EREL,ECOR,EFIN,DS,HMIN,HMAX,
     + MXSTEP,MXDG,MPREF,NFAC,NEWTON,ITOK,TOLINT,MXDM,D,IPIVOT,WCP,WP,
     + TANGNT,DIR,U,V,STO,STO2,PHI,PHIA,ISING,IRCODE)
      PARAMETER(MAXDEG=13,TOLIT=1.5,TOLF=1.5,NFCOR=15)
C***********************************************************************
C *** A B C O N : PREDICTOR-CORRECTOR CONTINUATION USING ADAMS-BASHFORTH
C ***             VARIABLE ORDER PREDICTORS  WITH ERROR-STEPSIZE CONTROL
C ***              (SINGLE PRECISION VERSION: UPDATED 5-27-93.)
C **********************************************************************
C The authors of this software are Bruce N. Lundberg and Aubrey B. Poore
C Permission to use, copy, modify, and distribute this software for any
C purpose without fee is hereby granted, provided that this entire notice
C is included in all copies of any software which is or includes a copy
C or modification of this software and in all copies of the supporting
C documentation for such software.
C THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
C WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
C REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
C OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
C
C E-MAIL:  lundberg@math.colostate.edu
C
C 1 ********  P U R P O S E  ******************************************
C
C THIS PROGRAM COMPUTES A SEQUENCE OF SOLUTIONS (UK,RK)=W(NP1)
C OF F(U,R) = 0  (WHERE F IS A FUNCTION FROM NP1=N+1 SPACE
C TO N SPACE) STARTING WITH A GIVEN INITIAL SOLUTION
C W = (UZERO,RZERO),AND TERMINATING AT RK = RLOW OR RHIGH,
C FINAL VALUES OF THE PARAMETER R.
C THE DETERMINANT OF THE SYSTEM DERIVATIVE IS MONITORED FOR
C CHANGES IN SIGN TO DETECT SINGULAR POINTS.
C
C THIS IS AN IMPLEMENTATION OF THE ARCLENGTH CONTINUATION
C ALGORITHM GIVEN IN THE PAPER :
C
C REFERENCE:  "VARIABLE ORDER ADAMS-BASHFORTH PREDICTORS WITH
C              ERROR-STEPSIZE CONTROL FOR CONTINUATION METHODS",
C              BY BRUCE N. LUNDBERG AND AUBREY B. POORE,
C              SIAM J. OF SCI. AND STAT. COMPUTATION (12)
C              P. 695-723, MAY 1991.
C
C 2 ******* CALL ARGUMENTS FOR ABCON **********************************
C
C ** INITIAL INPUT VARIABLES, FUNCTIONS, AND CONTROL PARAMETERS **
C N -- THE SIZ OF THE PROBLEM = NUMBER OF COMPONENTS IN F, THE 
C      CONTINUATION PROCEEDS IN THE N+1 VARIABLES W = (Z, R)
C WC(N+1) -- ON INPUT, CONTAINS INITIAL SOLUTION POINT TO F(W) = 0.
C       ON OUTPUT, CONTAINS LAST CORRECTED SOLUTION POINT.
C FFIN(N) -- ON OUTPUT, CONTAINS F AT FINAL POINT WC 
C F, DZF, DRF -- USER WRITTEN SUBROUTINES EVALUATION THE PROBLEM
C          FUNCTIONS AND DERIVATIVES.  NAMES MUST BE DECLARED EXTERNAL
C          IN THE CALLING PROGRAM.  (SEE SECTION 6 BELOW).
C NPRNT -- PRINT CONTROL CODE:
C    NPRNT=0, NO PRINTING OF INFORMATION
C    NPRNT=-3, PERFORMANCE COUNTS ONLY
C    NPRNT=-1, PRINT CORRECTED POINTS ONLY.
C    NPRNT=1, INITIAL HEADING WITH PARAMETER SETTINGS PRINTED
C             FIRST TWO STEPS PRINTED ALONG WITH SINGULAR POINT
C             INFO AND FINAL POINT INFO
C    NPRNT=2, ABOVE PLUS CORRECTED POINTS (FIRST TWO, EVERY KSKIP STEP)
C    NPRNT=3, ABOVE PLUS HEADING FOR STEPS AND OTHER STEP INFO.
C    NPRNT=4, ABOVE PLUS STATUS MESSAGES
C    NPRNT>4, ABOVE PLUS FINAL JACOBIAN CALCULATED AND PRINTED
C  K -- THE STEP NUMBER
C  DS -- THE SIZE OF THE PREDICTOR STEP TAKEN IN STEP K (K=1,2,..)
C  WC -- THE CORRECTED K TH TERM (UK,RK) OF THE GENERATED SEQUENCE
C           OF SOLUTIONS.
C  M+1 -- PREDICTOR USED IS AN A-B (M+1) STEP METHOD.
C  KNTROL -- THE NUMBER OF CORRECTOR ITERATIONS USED IN CORRECTION.
C
C KSKIP -- PRINT SKIPPING PARAMETER: PRINT EVERY KSKIP STEP.
C ORIENT -- SETS ORIENTATION FOR INITIAL DIRECTION OF CONTINUATION
C    ORIENT = +1, THE PARAMETER WC(N+1) INITIALLY IS INCREASING
C     ORIENT = -1, THE PARAMETER INITIALLY DECREASES
C NCOMP -- SETS THE TARGET COMPONENT FOR STOPPING THE CONTINUATION:
C    NCOMP = 1, ..., N+1, THE CONTINUATION WILL STOP AND LOCATE A POINT
C            WITH WC(NCOMP) = VCOMP ONCE WC(NCOMP) .GT. VCOMP
C    NCOMP = -1, ...,-(N+1) THE CONTINUATION WILL STOP AND LOCATE A
C            POINT WITH WC(-NCOMP) = VCOMP ONCE WC(-NCOMP) .LT. VCOMP
C    NCOMP = 0, THIS FEATURE IS DISABLED (SEE RLOW, RHIGH HOWEVER)
C VCOMP -- TARGET VALUE FOR WC(ABS(NCOMP)) (SEE NCOMP)
C RHIGH, RLOW -- THE BOUNDS OF CONTINUATION: IF THE PREDICTED OR
C                CORRECTED POINT CROSSES THE HYPERPLANES R = RHIGH OR
C                R = RLOW, THE CONTINUATION WILL STOP WITH A POINT WC
C                WITH WC(N+1) = RHIGH, (OR RLOW).
C KTARG --  MINIMUM NUMBER OF STEPS BEFORE FINAL POINT STOPPING WILL BE
C           INVOKED.  THIS CAN BE USED TO PICK OUT SUBSEQUENT PASSINGS
C           OF A PARAMETER VALUE FOR ATTENTION.  
C IFIN -- CONTROL PARAMETER FOR FINAL POINT CORRECTION
C    IFIN .GT. 0,  FINAL CORRECTION DONE WITH A NEWTON OR CHORD
C         ITERATION WITH TARGET PARAMETER FIXED, NO BORDERING
C         AND .LE IFIN ITERATIONS.  
C    IFIN .LT. 0, USES ROUTINE COREX BUT WITH DIR RESET TO GIVE
C         CORRECTION WITH TARGET PARAM. FIXED AND .LE. -IFIN ITERATIONS
C    IFIN .EQ 0,  USES ROUTINE COREX; BORDERING ALGORITHM
C              WITH ITERATES IN A HYP. PLANE ORTH. TO USUAL PREDICTOR.
C              MAXIMUM ITERATES DETERMINED AS WITH POINTS (SEE NEWTON)
C EABS, EREL -- RELATIVE AND ABSOLUTE PREDICTOR ERROR PARAMETERS
C               (DEFAULT VALUES COMPUTED FOR EACH ONE .LE. 0 ON ENTRY)
C ECOR -- ERROR PARAMETER FOR CORRECTED POINTS
C               (DEFAULT VALUE COMPUTED IF .LE. 0 ON ENTRY)
C EFIN -- ERROR PARAMETER FOR FINAL TARGET POINT:
C     EFIN .GT. 0, USE A CHORD CORRECTOR, STOPPING ERROR EFIN
C     EFIN .LT. 0, USE A NEWTON CORRECTOR, STOPPING ERROR -EFIN
C     EFIN .EQ. 0, USE NO CORRECTOR
C DS -- INITIAL STEPSIZE FOR FIRST PREDICTOR STEP
C               (DEFAULT VALUE COMPUTED IF .LE. 0 ON ENTRY)
C HMIN, HMAX -- MINIMUM AND MAXIMUM PREDICTOR STEPSIZES
C MXSTEP -- MAXIMUM NUMBER OF CONTINUATION STEPS (POINTS)
C MXDG -- THE MAXIMUM DEGREE OF THE INTERPOLATING POLYNOMIAL ALLOWED
C         IN THE PREDICTION PHASE (MXDG + 1 STEP METHOD).  THIS MAY
C         BE SPECIFIED BY THE USER FROM 0 TO 12, OR 
C     < 0, DEFAULT MXDG COMPUTED AS FOLLOWS:
C     = IMIN(12,IMAX(2,NINT(-ALOG10((EREL+EABS)/2.E0)))
C     THAT IS THE CLOSEST INTEGER TO THE NEGATIVE OF THE
C     COMMON LOG OF THE AVERAGE OF THE PREDICTOR ERRORS,
C     BUT AT MOST 12 AND AT LEAST 2.
C MPREF -- AUTOMATIC PREDICTOR ORDER INCREASE UNTIL ORDER > MPREF
C NFAC -- CONTROL PARAMETER FOR TANGENT COMPUTATION:
C     NFAC = 0 FOR COMPUTATION OF THE TANGENT AT WP (PREDICTED POINT)
C     NFAC = 1 FOR COMPUTATION OF THE TANGENT AT WC (CORRECTED POINT)
C NEWTON -- CONTROL PARAMETER FOR CORRECTOR:
C     NEWTON = 1, 2, ...  CORRECTOR ABORTS AFTER NEWTON CHORD ITERATES
C     NEWTON = -1, -2, ... CORRECTOR USES UP TO 'NEWTON' NEWTON ITERATES
C                          (JACOBIAN REEVALUATED FOR EACH ITERATE)
C ITOK --  MAX CORRECTOR ITERATES WHEN INCREASING NEXT STEPSIZES AFTER
C          CORR. FAILS(ITOK .GE.ABS(NEWTON) DISABLES THIS FEATURE).
C TOLINT -- ERROR TOLERANCE FOR MAX NORM OF F(W) AT INITIAL POINT
C           IF THIS TOLERANCE IS NOT MET, IRCODE = -6 RETURNS.
C
C  ** STORAGE AND WORK SPACE ARGUMENTS **
C MXDM -- LEADING DIMENSION OF CERTAIN TWO DOUBLE INDEXED WORK ARRAYS
C
C D(MXDM,N+1) -- ARRAY FOR STORING THE NXN JACOBIAN OF F WRT Z 
C IPIVOT(N+1) -- STORES PIVOT INFORMATION FOR LINEAR SYSTEM SOLVER
C WCP(N+1) -- STORAGE FOR PENULTIMATE CORRECTED POINT
C WP(N+1) -- PREDICTED SOLUTION POINT
C TANGNT(N+1) -- UNIT TANGENT TO PATH
C DIR(N+1) -- PREDICTION DIRECTION USED IN COREX AND PREDIC
C U(N+1), V(N+1), STO(N+1), STO2(N+1) -- STORAGE VECTORS FOR COREX
C PHI(MXDM,0:MAXDEG), PHIA(MXDM,0:MAXDEG) SCALED DIVIDED DIFF. STORAGE
C
C ISING -- INTEGER FLAG WHICH TURNS ON BIFURCAION POINT CHECK
C       = 0, NO CHECKING FOR BIFURCATION, ALL SINGULARITIES ASSUMED 
C            TO BE FOLDS
C       > 0, ROUTINE SINGCK USED TO ANALYSE SINGULAR POINTS
C ** ERROR RETURN CODE **
C IRCODE -- RETURNS INTEGER ERROR CODE:
C       =  0, NO ERRORS, NORMAL RETURN
C       = -1, MAXIMUM NUMBER OF STEPS COMPUTED WITHOUT REACHING TARGET
C       = -2, FAILURE IN FINAL PREDICTION
C       = -3, FAILURE IN FINAL CORRECTION
C       = -4, SINGULAR JACOBIAN OR ZERO PIVOT ENCOUNTER
C       = -5, INPUT ERROR, HMAX < HMIN
C       = -6, INPUT ERROR, INITIAL POINT DOES NOT SATISFY F=0 TO TOLINT
C       = -7, ZERO PIVOT IN INITIAL PLU
C       = K>0, CORRECTOR FAILED NFCOR TIMES AT STEP K
C
C 3 ********* INTERNAL FORTRAN PARAMETERS ******************************
C  (THESE ARE SET INTERNALLY BUT MAY BE CHANGED TO SAVE STORAGE, ETC.)
C MAXDEG --  THE MAXIMUM DEGREE OF INTERPOLATING POLYNOMIAL IN
C            THE PREDICTOR. AN M TH DEGREE POLY GIVES AN ADAMS-
C            BASHFORTH M+1 STEP METHOD. THE DEFAULT IS 13
C TOLIT, TOLF -- MINIMUM CONTRACTION PARAMETERS FOR COREX, SEE COREX
C NFCOR  --  THE MAXIMUM NUMBER OF CORRECTOR FAILURES ALLOWED
C            PER STEP; IF A CORRECTOR FAILS MORE THAN NFCOR
C            TIMES IN A STEP THE RUN IS ABORTED.
C
C 4 ********** LOCAL VARIABLES AND CONSTANTS ***************************
C NP1 -- = N + 1 AND EQUALS THE DIMENSION OF (U,R) IN F(U,R) = 0.
C M -- CURRENT DEGREE OF INTERP POLY(M+1 AND M+2 STEP METHODS)
C
C ** LOGICAL VARIABLES **
C    ORFLG -- =, .TRUE. IF DEGREE REDUCTION IS INDICATED IN SELST2.
C             =, .FALSE. IF DEGREE REDUCTION IS NOT INDICATED.
C     MXFLG -- =, .TRUE. IF CONVERGENCE IS NOT OBTAINED IN SELST2.
C              =, .FALSE. IF CONVERGENCE IS OBTAINED.
C
C 5 ********** FUNCTIONS TO BE WRITTEN BY THE USER *********************
C    THESE MUST BE CODED AS SUBROUTINES AND DECLARED EXTERNAL IN THE
C    CALLING PROGRAM AND HAVE THERE NAMES INCLUDED IN PLACE OF THE
C    ARGUMENTS F, DZF, DRF
C
C F -- THE FUNCTION TO BE CONTINUED FOR. (F(N,VALUE,ARGUMENT))
C DZF -- THE DERIVATIVE(NXN MATRIX) OF F WITH RESPECT TO THE FIRST
C        N VARIABLES.(DZF(N,MXDM,RESULT,ARGUMENT))
C DRF -- THE DERIVATIVE(N - VECTOR) OF F WITH RESPECT TO THE LAST
C        VARIABLE, THE PARAMETER R.(DRF(N,RESULT,ARGUMENT))
C
C 6 ********** SUBROUTINES USED ****************************************
C ** PREDICTOR SUBROUTINES **
C   PREDIC -- COMPUTES PREDICTOR , UPDATES, SELECTS ORDER AND STEPSIZE.
C              PREDICTOR POLYNOMIAL.
C       PSDIR -- PREDICTS SOLUTION AND PRED. DIRECTION.
C       PRDINF -- COMPUTES PREDICTOR INFORMATION.
C       SELST2 --- SELECTS STEPSIZE TO SPECIFIED TOLERANCE.
C       ARCLG2 -- COMPUTES AN ESTIMATE OF ARCLENGTH FOR UPDATING DS.
C       PDNRM --  USED IN ARCLG2 TO ESTIMATE TANGENT NORM.
C ** CORRECTOR SUBROUTINES **
C   COREX -- TO CORRECT THE RESULTS OF PREDICTION, APPROX. SOLUTION TO
C            F = 0.  ALSO COMPUTES NEW DW/DS AND SYSTEM DETERMINANT.
C       SGEFA,SGESL,SGEDI, --    SQUARE LINEAR SYSTEM SOLVER FROM
C                                 LINPACK..
C       DWDS -- COMPUTES THE UNIT TANGENT TO THE PATH AT REGULAR POINTS.
C       FINALP -- CHECKS FOR PASSAGE OF A TARGET POINT.  PREDICTS AND 
C               CORRECTS (OR SETS UP FOR A FINAL CALL OF COREX)
C               FINAL TARGET POINT.
C ** SINGULARITY CHECKING AND ANALYSIS **
C   SINGCK -- CHECKS FOR THE PASSAGE OF FOLD AND BIFURCATION POINTS.
C       MAKMAT -- USED IN SINGCK TO CONSTRUCT SYSTEM MATRIX.
C ** LINPACK ROUTINES USED **
C   SDOT   -- SINGLE PRECISION DOT PRODUCT.
C   SNRM2  -- SINGLE PRECISION EUCLIDEAN NORM.
C   SGEFA  -- SINGLE PRECISION PLU DECOMPOSITION
C   SGESL  -- SINGLE PRECISION SYSTEM SOLVER (USING PLU)
C   SGEDI  -- SINGLE PRECISION DETERMINANT COMPUTATION
C   SAXPY,SASUM, SSWAP, SSCAL  -- ROUTINES CALLED BY ABOVE
C 7 ****************** TYPE AND DIMENSION ***************************
      REAL WC(N+1),WCP(N+1),WP(N+1),DIR(N+1),U(N+1),FFIN(N+1),
     +    V(N+1),STO(N+1),STO2(N+1),D(MXDM,N+1),TANGNT(N+1),
C ARRAYS FOR PREDICTOR AND SCALED DIVIDED DIFFERENCE STORAGE
     +     G(0:MAXDEG,1:MAXDEG),PSYP(0:MAXDEG),PSY(0:MAXDEG),
     +     A(0:MAXDEG),B(0:MAXDEG),PHI(MXDM,0:MAXDEG),
     +     PSYARC(0:MAXDEG),PHIA(MXDM,0:MAXDEG),
C REAL SPECIFICATION VARIABLES
     +     EREL,EABS,ECOR,EFIN,HMIN,HMAX,RLOW,RHIGH,ORIENT,
C REAL STORAGE VARIABLES
     +     DS,DSP,DSM,DSMM1,DSMM2,DETOLD,DETNEW,TEMP
C
      INTEGER IPIVOT(N+1),IFLAG,KNTROL,LEX,MXDG,MXSTEP,IRCODE,
     +        NPRNT,M,N,NEWTON,NFAC,IFIN,NP1,MPREF,COUNT,ITOK,ISING
C
      LOGICAL ORFLG,FAIL,UPDATE,SELECT,FINAL
      COMMON KJAC,KFUNC,KLIN
      EXTERNAL F,DZF,DRF
C
C 8 *********** EXECUTABLE CODE BEGINS HERE **************************
C
C ****************** INITIALIZATION **********************************
C
C *** INITIALIZING  NP1,A,B,G ***
C
      IRCODE = 0
      NP1 = N + 1
      LEX = 0
      K = 1
      A(0) = 1.0
      B(0) = 1.0
      G(0,1) = 1.0
C ** TESTING ERROR PARAMETERS AND/OR COMPUTING DEFAULTS **
      EMACH = R1MACH(4)
C ** TESTING AND/OR COMPUTING DEFAULT STEPSIZE FOR FIRST STEP **
      IF(ECOR.LE.0.D0) ECOR = SQRT(EMACH)
      IF(EREL.LE.0.D0) EREL = EMACH**0.25D0
      IF(EABS.LE.0.D0) EABS = EMACH**0.25D0
      IF(ECOR.LE.1.E1*EMACH) ECOR = 1.E1*EMACH
      IF((EFIN.GT.0.D0).AND.(EFIN.LE.1.E1*EMACH)) EFIN = 1.E1*EMACH
      IF((EFIN.LT.0.D0).AND.(EFIN.GT.-1.E1*EMACH)) EFIN = -1.E1*EMACH
      IF(EREL.LE.1.E1*EMACH) EREL = 1.E1*EMACH
      EA = EABS
      ER = EREL
      EC = ECOR
C ** COMPUTING THE DEFAULT INITIAL STEPSIZE DS IF DS.LE.0 **
      IF(DS.LE.0.D0) THEN
         TEMP = SNRM2(NP1,WC,1)
         DS = 0.5*(SQRT(EABS) +SQRT(EREL)*TEMP )
      END IF
C ** COMPUTING THE DEFAULT MAXIMUM DEGREE (FOR PREDICTOR)
      IF(MXDG.LT.0) THEN
        MXDG = 1+NINT(-LOG10(.5E0*(EREL+EABS)))
        MXDG = MAX(MXDG,2)
      END IF
        MXDG = MIN(MXDG,12)
C  COUNTERS INITIALIZED
      KFUNC = 0
      KJAC = 0
      KLIN = 0
      KCUT = 0
      KSING = 0
      MPREF = MAX(MPREF,1)
      FAIL = .FALSE.
      FINAL = .FALSE.
C
C *** ECHO PRINT OF SPECIFICATIONS ************************************
C
C IF NPRNT .LT.1 SKIP PRINT
      IF(NPRNT.GT.0) THEN
      PRINT*,'********************************************************'
      PRINT*,'****** ARCLENGTH CONTINUATION ALONG A SOLUTION *********'
      PRINT*,'******        CURVE OF  F(U,R) = 0.            *********'
      PRINT*,'********************************************************'
      PRINT*,'CONTINUATION IN THE INTERVAL (',RLOW,' , ',RHIGH,')'
      PRINT*,'WITH INITIAL ORIENTATION: ',ORIENT
      PRINT*,'AND MAXIMUM NUMBER OF STEPS: ',MXSTEP
      PRINT*,' PRED STEPSIZE BOUNDS: MIN = ',HMIN,' MAX = ',HMAX
      PRINT*,' INITIAL STEPSIZE: ',DS
      PRINT*,' ERROR REQUIREMENTS FOR PREDICTOR: REL,ABS',EREL,EABS
      PRINT*,'THE MAX STEP METH. OF PRED: ',MXDG+1,'MPREF=',MPREF
      PRINT*,'PREFERED MAX. # CORR ITERATIONS ITOK=  ',ITOK
      PRINT*,'ERROR REQ. FOR CORR.: STEPS',ECOR,'FINAL POINT',ABS(EFIN)
      IF(NEWTON.LT.0) THEN
        PRINT*,'CORRECTOR USES ** NEWTONS METHOD **'
      ELSE
        PRINT*,'CORRECTOR USES ** CHORD METHOD **'
      END IF
      IF(NFAC.EQ.1) THEN
        PRINT*,'WITH DW/DS EVALUATED AT THE ** CORRECTED ** PT.'
      ELSE
        PRINT*,'WITH DW/DS EVALUATED AT THE ** PREDICTED ** PT.'
      END IF
      IF(NCOMP.NE.0) THEN
        PRINT*,'TARGET POINT: POINT WITH',ABS(NCOMP),'COMPONENT',VCOMP
      END IF
      PRINT*,'FINAL(TARGET) POINT COMPUTED WITH:'
      IF(NCOMP.GT.0) THEN
        IF(EFIN.LT.0.) PRINT*,'NEWTON ITERATION'
        IF(EFIN.GT.0.) PRINT*,'CHORD ITERATION'
        PRINT*,'WITH THE TARGET PARAMETER HELD FIXED'
      END IF
      IF(NCOMP.LT.0) THEN
        IF(EFIN.LT.0.) PRINT*,'NEWTON ITERATION'
        IF(EFIN.GT.0.) PRINT*,'CHORD ITERATION'
        PRINT*,'WITH BORDERING ALGORITHM AND TARGET PARAMETER FIXED'
      END IF
      IF(NCOMP.EQ.0) THEN
        IF(EFIN.LT.0.) PRINT*,'NEWTON ITERATION'
        IF(EFIN.GT.0.) PRINT*,'CHORD ITERATION'
        PRINT*,'WITH ITERATES IN A HYPERPLANE ORTH. PRED. DIRECTION'
      END IF
      IF(EFIN.EQ.0) THEN
        PRINT*,'WITH NO CORRECTION (PREDICTION RETURNED)'
      END IF
       PRINT*,'############### START OUTPUT #####################'
      END IF
C
C **************** TEST DATA ************************************
C
      IF(HMIN.GT.HMAX) THEN
       IF(NPRNT.GT.3) THEN
        PRINT*,' ** DATA ERROR: STEPSIZE BOUNDS **'
       END IF
       IRCODE = -5
                                            RETURN
      END IF
      CALL F(N,STO,WC)
      TEMP = 0.0
      DO 25 L=1,N
         TEMP = MAX(TEMP,ABS(STO(L)))
 25   CONTINUE
      IF(TEMP.GT.TOLINT) THEN
       IF(NPRNT.GT.1) THEN
        PRINT*,'DATA ERROR: INITIAL POINT NOT A ZERO OF F(X,R)'
        PRINT*,'CHANGE ERROR TOLERANCE IN CALLING PROGRAM'
        PRINT*,'OR OBTAIN MORE ACCURACY IN INITIAL POINT'
        PRINT*,'F : ',(STO(L),L=1,N)
        PRINT*,'SUP-NORM(F): ',TEMP
       END IF
       IRCODE = -6
                                            RETURN
      END IF
C
C ***************** INITIAL STEP *************************************
C
      M = 0
      IFAILF = 0
      ARCLTH = 0.
      ARCP = 0.
      ARCS = 0.
      CALL DZF(N,MXDM,D,WC)
      CALL DRF(N,STO,WC)
      CALL SGEFA(D,MXDM,N,IPIVOT,IFLAG)
      IF(IFLAG.GT.0) THEN
        IF(NPRNT.GE.3) THEN
          PRINT*,'ZERO PIVOT IN INITIAL PLU'
        END IF
        IRCODE = -7
                                       RETURN
      END IF
      CALL SGEDI(D,MXDM,N,IPIVOT,DETOLD,V,10)
      CALL DWDS(D,N,MXDM,STO,TANGNT,IPIVOT,ORIENT)
28    DO 30 L=1,NP1
         PHI(L,0) = TANGNT(L)
         WCP(L) = WC(L)
30    CONTINUE
C
C** FIRST PREDICTION-CORRECTION **
C
      COUNT = 0
      CALL PSDIR(M,N,MAXDEG,MXDM,DS,WCP,WP,DIR,G,B,PHI)
C
C  CHECK BOUNDS OF CONTINUATION AND REDUCE STEPSIZE SO FIRST STEP
C  IS NOT THE LAST STEP
C
      IF(WP(NP1).LE.RLOW) DS =0.5*(RLOW-WC(NP1))/TANGNT(NP1)
      IF(WP(NP1).GE.RHIGH) DS =0.5*(RHIGH-WC(NP1))/TANGNT(NP1)
      IF(NCOMP.GT.0) THEN
        IF((WP(NCOMP).GE.VCOMP).AND.(WC(NCOMP).LE.VCOMP)) THEN
          DS=0.5*(VCOMP-WC(NCOMP))/TANGNT(NCOMP)
        END IF
      ELSE IF(NCOMP.LT.0) THEN
        IF((WP(-NCOMP).LE.VCOMP).AND.(WC(-NCOMP).GE.VCOMP)) THEN
          DS=0.5*(VCOMP-WC(-NCOMP))/TANGNT(-NCOMP)
        END IF
      END IF
C
32    CALL PSDIR(M,N,MAXDEG,MXDM,DS,WCP,WP,DIR,G,B,PHI)
C
         KNTROL = NEWTON
         CALL COREX(N,MXDM,WC,WP,DIR,TANGNT,D,U,V,STO,STO2,EC,ORIENT,
     +    DETNEW,KNTROL,NFAC,IPIVOT,ENRM,TOLIT,TOLF,F,DZF,DRF,NPRNT)
C
C ** HALVING STEPSIZE IF CORRECTOR FAILS -- TRY NFCOR TIMES **
         IF(KNTROL.EQ.0)  THEN
           IF(NPRNT.GT.3) PRINT*,'** CORRECTOR FAILURE **'
          IF(COUNT.LT.NFCOR) THEN
            COUNT = COUNT + 1
            DS = 0.5*DS
            GO TO 32
          ELSE
            IRCODE = 1
                                             RETURN
          END IF
       END IF
C
C ** INITIAL OUTPUT **************************************************
C
         IF(ABS(NPRNT).EQ.1) THEN 
           PRINT*,(WCP(L),L=1,NP1)
           PRINT*,(WC(L),L=1,NP1)
          ELSE IF(NPRNT.GT.2) THEN
           PRINT*,'THE INITIAL POINT: (U,R): ',(WCP(L),L=1,NP1)
           PRINT*,'** STEP 1 **: (U,R): ',(WC(L),L=1,NP1)
51         PRINT*,'TANGENT: ',(TANGNT(L),L=1,NP1)
           PRINT*,'SUPNORM(WC-WP):',ENRM
         END IF
C
55       CONTINUE
C
         PSYP(0) = DS
         PSY(0) = DS
         KPRNT = KSKIP
C
C*********************************************************************
C ** MAIN LOOP: PREDICTOR-CORRECTOR STEPS **
C*********************************************************************
C
      DO 1000 K=2,MXSTEP
C ** SHIFTING  AND  UPDATING **
         DETOLD = DETNEW
         DSP = DS
         ARCP = ARCP + DS
         UPDATE = .TRUE.
         SELECT = .TRUE.
         COUNT = 0
C
C ******* P R E D I C T I O N ***************************************
C
         IF(FAIL) THEN
           DS = DSP
           SELECT = .FALSE.
        ELSE
          SELECT = .TRUE.
          EA = EABS
          ER = EREL
        END IF
C
200     CALL PREDIC(M,N,MAXDEG,MXDM,WCP,WC,TANGNT,DSP,DS,
     +   WP,DIR,PSY,PSYP,PSYARC,PHI,PHIA,G,A,B,STO,EA,ER,
     +   HMIN,HMAX,SELECT,UPDATE,LEX,MPREF,MXDG,K)
C
         IF(NPRNT.GT.1) THEN
           PRINT*,'PRED. STP: ',DSP,' CORR. ITER',KNTROL
           PRINT*,'************** ',K,' ****************************'
         END IF
        KNTROL = NEWTON
C  ** UPDATE ARCLENGTH ACCUMULATION **
         IF(UPDATE) ARCLTH = ARCLTH + DSP
C
C?????????  B O U N D   A N D   T A R G E T   C H E C K I N G ??????????
C
C ** CHECKING IF PREDICTOR IS OUTSIDE OF THE BOUNDS OF CONTINUATION **
C    SET UP (OR PREFORM) FINAL CORRECTION AS INDICATED (FINAL = .TRUE.)
        IF(K.GE.KTARG) THEN
          CALL FINALP(N,M,MXDM,MAXDEG,D,WCP,WP,V,U,DS,DIR,EFIN,
     +         RLOW,RHIGH,LEX,KNTROL,G,PSYP,PSY,A,B,PHI,IPIVOT,WC,
     +         NCOMP,VCOMP,F,DZF,DRF,NPRNT,FINAL,IFIN,EC,IRCODE)
          IF(FINAL) THEN
           IF((KNTROL.EQ.0).AND.(EFIN.NE.0.))  GO TO 225
           IF(EFIN.EQ.0.) THEN
            IF(NPRNT.GT.3) PRINT*,'NO FINAL CORRECTION'
            KNTROL = 0
            GO TO 1890
           END IF
           IF((IFIN.GT.0).AND.(KNTROL.GT.0)) THEN
            IF(NPRNT.GT.3) PRINT*,'CORRECTION IN FINALP'
            GO TO 889
           END IF
           IF((IFIN.LT.0).AND.(NPRNT.GT.3)) THEN
            PRINT*,'VERTICAL CORRECTION IN COREX'
           END IF
           IF((IFIN.EQ.0).AND.(NPRNT.GT.3)) THEN
            PRINT*,'HYPERPLANE CORRECTION, FINAL'
           END IF
          END IF
        END IF
C
C *******  C O R R E C T I O N  *************************************
C
        CALL COREX(N,MXDM,WC,WP,DIR,TANGNT,D,U,V,STO,STO2,EC,ORIENT,
     +   DETNEW,KNTROL,NFAC,IPIVOT,ENRM,TOLIT,TOLF,F,DZF,DRF,NPRNT)
C
C ??????? CORRECTOR FAILURE ??????????????????????
C
C ** THE FOLLOWING CODE HANDLES CORRECTOR FAILURES
C     REDUCING STEPSIZE IF CORRECTOR FAILS -- TRY NFCOR TIMES **
C
225         CONTINUE
            IF(KNTROL.EQ.0)  THEN
              FINAL = .FALSE.
              EC = ECOR
              IF(NPRNT.GT.3) PRINT*,'   ** CORRECTOR FAILURE **'
              KCUT = KCUT + 1
              SELECT = .FALSE.
              UPDATE = .FALSE.
              IF(COUNT.LT.NFCOR) THEN
               COUNT = COUNT + 1
               IF((COUNT.EQ.1).AND.(.NOT.FAIL)) THEN
                M = 0
                EA = .01*EABS
                ER = .01*EREL
                CALL SELST2(M,N,MAXDEG,MXDM,DS,DSP,EA,ER,HMIN,HMAX,
     +           WCP,WP,G,PSYP,PSY,PHI,B,A,DIR,ORFLG)
                DS = MIN(DS,DSP)
               ELSE
                DS = 0.5*DS
               END IF
               GO TO 200
              ELSE
               IRCODE = K
                                              RETURN
              END IF
            END IF
C ** DECIDING IF NEXT STEP SHOULD BE CAUTIOUS  RESET COUNT **
         IF(KNTROL.LE.ITOK) FAIL = .FALSE.
         IF(COUNT.GT.0) FAIL = .TRUE.
C
C ??????????? SINGULARITY DETECTION  ????????????????????????????
C
255      IF(DETNEW*DETOLD.LE.0.0) THEN
           IF(ISING.GT.0) THEN
           CALL SINGCK(N,MXDM,WC,WCP,D,PHI(1,0),PHI(1,0),V,DZF,DRF,
     +     INDIC,IPIVOT)
           ELSE
           INDIC = 1
           END IF
           IF(INDIC.GT.0) THEN
             KSING = KSING + 1
             IF(NPRNT.GT.1) THEN
               PRINT*,'ARCLENGTH FROM LAST SING: ',ARCLTH-ARCS
             END IF
             ARCS = ARCLTH
           END IF
            IF(INDIC.EQ.1) THEN
             IF(NPRNT.GT.0) THEN
             PRINT*,'*** A FOLD POINT HAS BEEN PASSED ***'
            END IF
C ** CHANGE ORIENTATION AND REVERSE TANGENT **
              ORIENT = -ORIENT
              DO 275 L=1,NP1
                TANGNT(L) = -TANGNT(L)
275           CONTINUE
            END IF
         END IF
         IF((NPRNT.GT.0).and.(INDIC.EQ.2)) THEN
             PRINT*,'*** A BIFURCATION POINT HAS BEEN PASSED ***'
         END IF
         INDIC = 0
C
C ****  PRINTED OUTPUT ***********************************************
C
889      CONTINUE
         IF(FINAL.AND.(NPRNT.GT.1)) THEN
            PRINT*,'*******  FINAL CORRECTION  ********'
            KPRNT = KSKIP
         END IF
         IF((KPRNT.LT.KSKIP).OR.(KSKIP.LT.0)) THEN
           KPRNT = KPRNT+1
           GO TO 900
         ELSE
           KPRNT = 0
         END IF
C SKIP THIS PRINTING IF KPRNT<KSKIP ETC.
         IF(NPRNT.GT.1) THEN
           PRINT*,'(U,R) CORR.: = ',(WC(L),L=1,NP1)
         END IF
         IF(NPRNT.GT.2) THEN
           PRINT*,'PRED STP:',DS,M+1,'STP.METH. CORR.ITER:',KNTROL
           PRINT*,'WP =',(WP(L),L=1,NP1)
           IF((.NOT.FINAL).OR.(FINAL.AND.(IFIN.LE.0))) THEN
             PRINT*,'TANGENT:',(TANGNT(L),L=1,NP1)
             PRINT*,'SUPNORM(WC-WP):',ENRM
           END IF
         END IF
C
         IF(NPRNT.EQ.1) THEN
           PRINT*,(WC(L),L=1,NP1)
         END IF
900      CONTINUE
C
         IF(FINAL) GO TO 1890
1000  CONTINUE
C *****  E N D  OF  M A I N   L O O P ********************************
C
      IF(NPRNT.GT.3) THEN
      PRINT*,'**** MAX NUMBER OF STEPS TAKEN ****'
      END IF
      IRCODE = -1
      GO TO 5000
C
C ******** F I N A L   O U T P U T ***********************************
C
1890  CONTINUE
      CALL F(N,FFIN,WC)
      IF(NPRNT.GT.2) THEN
      TEMP = 0.0
      DO 1925 L=1,N
         TEMP = MAX(TEMP,ABS(FFIN(L)))
1925  CONTINUE
        PRINT*,'F : ',(FFIN(L),L=1,N)
        PRINT*,'SUP-NORM(F): ',TEMP
        PRINT*,'CORR.ITER.:  ',KNTROL
        PRINT*,'===================================================='
      END IF
5000  CONTINUE
      IF(ABS(NPRNT).GT.2) THEN
        PRINT*,'#PTS,JACBNS,FUNCTS,LIN SYS,ARCLGTH,CORTOL,SING,FAIL'
        PRINT*,K,KJAC,KFUNC,KLIN,ARCLTH+DS,ECOR,KSING,KCUT
      END IF
C PRINT JACOBIAN IF NPRNT .GT.3
            IF(NPRNT.GT.4) THEN
            CALL DZF(N,MXDM,D,WC)
            CALL DRF(N,STO,WC)
            PRINT*,'**********  FINAL JACOBIAN *********'
              DO 5278 LL=1,N
                 PRINT*,(D(LL,II),II=1,N),STO(LL)
5278           CONTINUE
            PRINT*,'********** END JACOBIAN *********'
            END IF
                                    RETURN
      END
C
C ********  P R E D I C  *********************************************
C
      SUBROUTINE PREDIC(MDEG,N,MAXDEG,MXDM,WCP,WC,TANGNT,DSP,DS,
     + WP,DIR,PSY,PSYP,PSYARC,PHI,PHIA,G,A,B,STO,EABS,EREL,HMIN,HMAX,
     + SELECT,UPDATE,LEX,MPREF,MXDG,K)
C*********************************************************************
C *** PURPOSE ***
C  THIS SUBROUTINE COMPUTES A PREDICTED SOLUTION WP(AND
C  AN ASSOCIATED PREDICTION DIRECTION DIR) TO A FIRST ORDER
C  SYSTEM OF O.D.E.'S, USING AN ADAMS-BASHFORTH PREDICTOR
C  SCHEME WITH A STEPSIZE AND ORDER CONTROL. 
C  PREDIC IS WRITTEN BY B. N. LUNDBERG FOR USE WITH THE SOFTWARE ABCON,
C  A FORTRAN SUBROUTINE FOR CONTINUATION ALONG A CURVE DESCRIBED IN:
C  "VARIABLE ORDER ADAMS-BASHFORTH PREDICTORS WITH AN ERROR-
C  STEPSIZE CONTROL FOR CONTINUATION METHODS",  BY B. N. LUNDBERG
C  AND A. B. POORE, SISSC (12), MAY 1991, P.695-723.
C  DATE: 5/26/89. 
C  THE  CODE PROCEEDS BY FIRST ESTIMATING ARCLENGTH,
C  UPDATING SCALED DIVIDED DIFFERENCES ETC.FROM THE
C  PREVIOUS STEP. THE ROUTINE PROVIDES THE OPTIONS
C  OF SKIPPING THE UPDATING  AND/OR THE STEP AND ORDER
C  SELECTION  OPTIONS OF THE CODE.
C  (SEE ALSO DOC. PROVIDED WITH THE ANCILARY
C   SUBROUTINES PSDIR,PRDINF,SELSTP.
C
C ** TYPE AND DIMENSION **
      REAL WC(N+1),WCP(N+1),WP(N+1),DIR(N+1),
     +     STO(N+1),TANGNT(N+1),
C ARRAYS FOR PREDICTOR AND SCALED DIVIDED DIFFERENCE STORAGE
     +     G(0:MAXDEG,1:MAXDEG),PSYP(0:MAXDEG),PSY(0:MAXDEG),
     +     A(0:MAXDEG),B(0:MAXDEG),PHI(MXDM,0:MAXDEG),
     +     PSYARC(0:MAXDEG),PHIA(MXDM,0:MAXDEG),
C REAL SPECIFICATION VARIABLES
     +     EREL,EABS,HMIN,HMAX,
C REAL STORAGE VARIABLES
     +     DS,DSP,DSM,DSMM1,DSMM2
C
C
      INTEGER MDEG,N,MAXDEG,MXDM,LEX,MPREF,MARC,K,MXDG
C
      LOGICAL ORFLG,SELECT,UPDATE
C
C*******************  EXPLANATION OF VARIABLES ***********************
C
C *** CONTROL VARIABLES ***
C
C  UPDATE,  = .TRUE. , IF UPDATE OF STEPSIZE,PSY,PSYP,PHI,WCP,DSP
C                      IS TO BE PERFORMED(I.E. ON THE FIRST CALL
C                      OF PREDICT IN A GIVEN STEP)
C             = .FALSE., OTHERWISE.
C
C  SELECT,  =  .TRUE. , IF STEPSIZE SELECTION(AND ORDER SELECTION)
C                       IS REQUESTED.
C           =  .FALSE., IF STEPSIZE AND ORDER SELECTION ARE NOT
C                       REQUESTED. IN THIS CASE THE DS AND MDEG
C                       GIVEN ON INPUT ARE USED IN PREDICTION.
C
C  MPREF --    ORDER SELECTION WILL BE PERFORMED ONLY AFTER MDEG
C              EXCEEDS THIS VALUE. MUST SET MPREF .GE.1
C
C  LEX ,    =  1 ,  IF LOCAL EXTRAPOLATION IS REQUESTED
C           =  0 ,  IF NO LOCAL EXTRAPOLATION IS REQUESTED.
C
C  K  -- THE STEP NUMBER(STARTING WITH K=1)
C
C  MXDG -- THE MAXIMUM VALUE OF MDEG(MUST BE SET .LT. MAXDEG)
C
C
C
C ** INPUT **
C   MDEG -- DEGREE OF INTERP POLY FOR WHICH STEPSIZE
C           SELECTION FOR PREDICTOR IS SOUGHT(MDEG.LE.MAXDEG-1)
C   N -- DIMENSION OF U IN THE PROBLEM F(U,R) = 0(N.LE.MXDM-1)
C   MXDM  -- THE MAX. DIMENSION OF U IN F(U,R), AS DIMENSIONED IN
C            THE CALLING PROGRAM.
C   MAXDEG -- MAXIMUM DEGREE INTERPOLATING POLY ALLOWED.
C   DS -- STEPSIZE TO BE USED(IN CASE SELECT = .FALSE)
C   EABS  -- ABSOLUTE ERROR REQUIREMENT FOR PREDICTOR
C   EREL  -- RELATIVE ERROR REQUIREMENT.
C   HMIN -- MINIMUM STEPSIZE ALLOWED.
C   HMAX -- MAXIMUM STEPSIZE ALLOWED.
C **** THE FOLLOWING ARE UPDATED ONLY ON THE
C   FIRST CALL PER STEP(WHEN UPDATE = .TRUE.)
C
C   DSP -- STEPSIZE USED IN PREVIOUS STEP.
C   WC  -- CORRECTED SOLUTION TO F = 0 FROM PREVIOUS STEP.
C   WCP  -- CORRECTED SOLUTION TO F = 0 FROM STEP BEFORE LAST
C   TANGNT  -- THE LAST UNIT TANGENT COMPUTED(AT WC OR PREVIOUS WP)
C   PSYP -- OLD VALUES OF PSY, FOR USE IN UPDATING.
C   DIR     -- THE PREDICTION DIRECTION COMPUTED IN THE PREVIOUS STEP
C         (FROM STEP BEFORE LAST)
C   B  -- VALUES OF B FROM LAST CALL OF PRDINF. FOR USE IN
C        ARCLENGTH UPDATING.
C   PSY  -- VALUES OF PSY FROM PREVIOUS STEP FOR USE
C           IN UPDATING.
C   PHI  -- SCALED DIVIDED DIFFERENCES USED IN CONSTRUCTING
C           PREDICTOR IN PREVIOUS STEP(FOR UPDATING TO NEW POLY)
C
C ** OUTPUT **
C   WP  --  THE PREDICTED SOLUTION POINT.
C   DIR --  THE PREDICTION DIRECTION: WP = WC + DS*DIR
C   MDEG  -- THE ORDER SELECTED
C   DSP -- THE PREVIOUS STEPSIZE IN ARCLENTH
C   PSY,PSYP,DSP,PHI,B,WCP  UPDATES.
C   DS -- THE STEPSIZE SELECTED(NOT CORRECTED TO ARCLENGTH)
C
C ** STORAGE AND ARRAYS USED BY OTHER SUBROUTINES CALLED HEREIN **
C   A,B,G,PSY,PSYP,PHI,PHIA,PSYARC  -- ARRAYS USED IN THE
C   COMPUTATION OF THE ADAMS-BASHFORTH PREDICTORS (THESE ARRAYS ARE USED
C   IN SUBROUTINE PRDINF AND PSDIR,SELSTP,ARCLG2,PDNRM)
C
C ** SUBROUTINES CALLED **
C    PSDIR -- PREDICTS SOLUTION AND PRED. DIRECTION.
C    PRDINF -- COMPUTES PREDICTOR INFORMATION
C    SELST2 --- SELECTS STEPSIZE TO SPECIFIED TOLERANCE.
C    ARCLG2 -- COMPUTES AN ESTIMATE OF ARCLENGTH FOR UPDATING
C              PREDICTOR POLYNOMIAL.
C    PDNRM -- USED IN ARCLG2.
C
C * SHIFTING  AND  UPDATING (IF UPDATE = .TRUE.) *
      IF(UPDATE) THEN
         IF(SELECT.AND.(K.GT.2).AND.(MDEG.LT.MXDG)) MDEG=MDEG+1
         IF(K.LE.1) THEN
           B(0)=1.0
           PSYP(0) = DSP
           DO 50 L=1,N+1
              PHI(L,0) = TANGNT(L)
50         CONTINUE
          GO TO 170
         END IF
        MARC = MAX(0,MDEG-1)
        IF(K.LE.2) THEN
          LEXP = 0
        ELSE
          LEXP = LEX
        END IF
        CALL ARCLG2(N+1,MXDM,MAXDEG,MARC,LEXP,DSP,TANGNT,DIR,PSYARC,
     +  PSY,PSYP,STO,PHI,PHIA,B,WC,WCP)
         DO 70 L=1,N+1
            WCP(L) = WC(L)
70       CONTINUE
         DO 80 J=0,MDEG
            PSYP(J) = PSY(J)
80       CONTINUE
C SHIFTING SCALED DIVIDED DIFFERENCES
         NN = MIN(MDEG+1,MAXDEG)
         DO 130 J=NN,1,-1
            DO 120 L=1,N+1
               PHI(L,J) = PHI(L,J-1)
120         CONTINUE
130      CONTINUE
         DO 140 L=1,N+1
            PHI(L,0) = TANGNT(L)
140      CONTINUE
C UPDATING SCALED DIVIDED DIFFERENCES
         DO 160 J=1,NN
            DO 150 L=1,N+1
               PHI(L,J) = PHI(L,J-1) - B(J-1)*PHI(L,J)
150         CONTINUE
160      CONTINUE
C
C
C
170   END IF
C ********** STEPSIZE AND ORDER SELECTION ****************
      LEX = 1
C
      IF(SELECT) THEN
            CALL SELST2(MDEG,N,MAXDEG,MXDM,DS,DSP,EABS,EREL,
     +      HMIN,HMAX,WCP,WP,G,PSYP,PSY,PHI,B,A,DIR,ORFLG)
CPRINT      PRINT*,'+++ DS SELECTED +++',DS
C
C **** ORDER SELECTION ********
C
C
          IF(MDEG.GE.MPREF) THEN
            DSM = DS
            CALL SELST2(MDEG-1,N,MAXDEG,MXDM,DSMM1,DSP,EABS,EREL,
     +      HMIN,HMAX,WCP,WP,G,PSYP,PSY,PHI,B,A,DIR,ORFLG)
CPRINT      PRINT*,'+++ DSMM1 SELECTED +++',DSMM1
C
        IF(MDEG.GE.2) THEN
            CALL SELST2(MDEG-2,N,MAXDEG,MXDM,DSMM2,DSP,EABS,EREL,
     +      HMIN,HMAX,WCP,WP,G,PSYP,PSY,PHI,B,A,DIR,ORFLG)
CPRINT      PRINT*,'+++ DSMM2 SELECTED +++',DSMM2
C
            IF((DSMM2.GE.DSM).AND.(DSMM2.GE.DSMM1)) THEN
              DS = DSMM2
              MDEG = MDEG-2
            ELSE
              IF(DSMM1.GE.DSM) THEN
                DS = DSMM1
                MDEG = MDEG-1
              END IF
            END IF
         ELSE
            IF((MDEG.EQ.1).AND.(DSMM1.GE.DSM)) THEN
              MDEG = 0
              DS = DSMM1
            END IF
         END IF
C
        END IF
C
       END IF
C ********** PREDICTION ******************************************
C
C  SELECTING LOCAL EXTRAP OR NOT
         IF(DS.LT.DSP) LEX = 0
C
C
200         CALL PRDINF(MDEG,MAXDEG,DS,PSYP,PSY,G,A,B)
C
         CALL PSDIR(MDEG+LEX,N,MAXDEG,MXDM,DS,WCP,WP,DIR,G,B,PHI)
C
C
C
                                                RETURN
1000  END
C
C
      SUBROUTINE PSDIR(MDEG,N,MAXDEG,MXDM,DS,WC,WP,DIR,G,B,PHI)
      REAL DS,WC(N+1),WP(N+1),G(0:MAXDEG,1:MAXDEG),
     +B(0:MAXDEG),PHI(1:MXDM,0:MAXDEG),DIR(N+1)
      INTEGER MDEG,N,MAXDEG,MXDM
C
C THIS SUBROUTINE COMPUTES THE PREDICTION DIRECTION AND SOLUTION
C POINT  ON A SOLUTION BRANCH OF F(U,R) = 0, USING AN ADAMS-
C BASHFORTH (MDEG + 1) STEP METHOD APPLIED TO THE DAVIDENKO O.D.E.
C WITH ARCLENGTH PARAMETERIZATION.
C PSDIR IS WRITTEN BY B. N. LUNDBERG FOR USE WITH THE SOFTWARE ABCON.
C    8/18/87
C ** INPUT **
C    DS -- THE STEPSIZE TO BE TAKEN
C    WC -- THE CURRENT CORRECTED SOLUTION TO F = 0
C    G,B -- COEFICIENTS FOR THE PREDICTOR COMPUTED BY A CALL
C           PRDINF
C    PHI -- SCALED DIVIDED DIFFERENCES
C    MXDM -- THE MAXIMUM DIMENSION OF THE PROBLEM: N.LE.MXDM-1.
C    MDEG -- DEGREE OF INTERPOLATING POLYNOMIAL USED IN PREDICTOR
C    MAXDEG -- MAXIMUM DEGREE OF INTERPOLATING POLYNOMIAL
C              ALLOWED(MAXDEG+1 IS THE MAXIMUM A-B STEP METHOD)
C              THUS MDEG .LE. MAXDEG-1.
C    N -- DIMENSION OF U IN THE PROBLEM F(U,R) = 0.(N.LE.MXDM-1)
C ** OUTPUT **
C    WP -- PREDICTED SOLUTION
C    DIR -- PREDICTION DIRECTION(NEEDED IN CORRECTION PHASE)
C CHECK MDEG
      IF(MDEG.GT.MAXDEG)  THEN
        PRINT*,'MDEG TOO LARGE IN PSDIR'
        STOP
      END IF
C
      DO 50 J=1,N+1
         DIR(J) = 0.0
50    CONTINUE
      DO 200 L=0,MDEG
         DO 100 J=1,N+1
            DIR(J) = DIR(J) + G(L,1)*B(L)*PHI(J,L)
100      CONTINUE
200   CONTINUE
C COMPUTING PREDICTOR
      DO 300 J=1,N+1
         WP(J) = WC(J) + DS*DIR(J)
300   CONTINUE
                              RETURN
      END
C
C
C
C ****************  P R D I N F  ************************************
C
C
      SUBROUTINE PRDINF(MDEG,MAXDEG,DS,PSYP,PSY,G,A,B)
      REAL DS,PSYP(0:MAXDEG),PSY(0:MAXDEG),
     +G(0:MAXDEG,1:MAXDEG),A(0:MAXDEG),B(0:MAXDEG)
      INTEGER MDEG,MAXDEG
C THIS SUBROUTINE COMPUTES THE INFORMATION NECESSARY FOR FORMING
C THE COEFFICIENTS USED IN COMPUTING THE (MDEG+1) AND
C (MDEG+2) STEP ADAMS-BASHFORTH METHODS.
C THIS INFORMATION IS NEEDED IN THE CALL TO SUBROUTINE PSDIR.
C SEE 'COMPUTER SOLUTION TO O.D.E.S' BY SHAMPINE AND GORDON.
C PRDINF IS WRITTEN BY B. N. LUNDBERG FOR USE WITH THE SOFTWARE ABCON.
C    8/18/87
C
C ** INPUT **
C    MDEG -- THE LOWER DEGREE OF THE INTERPOLATING POLYNOMIAL
C            TO BE USED IN THE ADAMS-BASHFORTH PREDICTOR:
C            (MDEG + 1)-STEP AND (MDEG+2)-STEP METHOD INFO IS
C             COMPUTED HEREIN. MUST HAVE MDEG.LE.MAXDEG-1 IN PRDINF.
C    MAXDEG -- THE MAXIMUM DEGREE INTERPOLATING POLY ALLOWED.
C    DS -- STEPSIZE USED IN PREDICTION.
C    PSYP -- PREVIOUS PSY VALUES.
C *** OUTPUT ***
C    G(0,J) -- INITIALIZED WITH 1/J,J=1,MAXDEG.
C    G(1,J) -- INITIALIZED WITH 1/J*(J+1),J=1,MAXDEG.
C    PSY -- NEW VALUES OF PSY CORRESPONDING TO STEPSIZE DS.
C           THESE MUST BE SHIFTED TO PSYP AFTER STEPSIZE SELECTION
C           IS COMPLETE.
C    A,B -- VALUES USED IN COMPUTING SCALED DIVIDED DIFFERENCES
C           AND A-B PREDICTOR RESPECTIVELY.
C    G(I,1) -- COEFFICIENTS FOR PREDICTOR(I=0,MDEG+1)
C ** STORAGE **
C    G -- REMAINING ENTRIES OF G USED AS WORKSPACE
C
C
C CHECKING MDEG.GT.MAXDEG.
      IF(MDEG.GE.MAXDEG) THEN
        PRINT*,'MDEG TOO LARGE IN PRDINF'
        STOP
      END IF
C INITIALIZING G(0-1,J)(G(0,MAXDEG+1)=1/(MAXDEG+1) BUT IS NOT USED)
      DO 50 J=1,MAXDEG
         G(0,J) = 1.0/FLOAT(J)
         G(1,J) = 1.0/FLOAT(J*(J+1))
50    CONTINUE
C COMPUTING PSY,B,A
      A(0) = 1.0
      B(0) = 1.0
      PSY(0) = DS
      DO 100 L=1,MDEG+1
         PSY(L) = PSYP(L-1) + DS
         B(L) = B(L-1)*PSY(L-1)/PSYP(L-1)
         A(L) = DS/PSY(L)
100   CONTINUE
C ** COMPUTATION OF G(I,1) **
C
      IF(MDEG.LT.1)            RETURN
      DO 300 I=2,MDEG+1
         DO 200 J=1,MDEG+2-I
            G(I,J) = G(I-1,J)-A(I-1)*G(I-1,J+1)
200      CONTINUE
300   CONTINUE
                                  RETURN
      END
C
C
C
C **************  S E L S T 2  ***************************************
C
C
      SUBROUTINE SELST2(MDEG,N,MAXDEG,MXDM,DS,DSP,EA,ER,HMIN,HMAX,
     +WC,WP,G,PSYP,PSY,PHI,B,A,DIR,ORFLG)
      PARAMETER(ALIMI=10.,ALIMD=.1)
      REAL DS,DSP,EA,ER,HMIN,HMAX,G(0:MAXDEG,1:MAXDEG),WC(N+1),
     +WP(N+1),PSYP(0:MAXDEG),PSY(0:MAXDEG),A(0:MAXDEG),
     +B(0:MAXDEG),DIR(N+1),PHI(MXDM,0:MAXDEG),
     +TEMP,WNRM,DIFNRM,ALPHA
      INTEGER MDEG,N,MAXDEG
      LOGICAL ORFLG
C
C ** PARAMETERS **
C    ALIMI --  LIMIT ON THE INCREASE IN STEPSIZE.
C    ALIMD --  LIMIT ON THE DECREASE IN STEPSIZE.
C
C THIS SUBROUTINE SELECTS A STEPSIZE TO ACHIEVE SPECIFIED ERROR
C REQUIREMENTS IN AN ADAMS-BASHFORTH MDEG+1STEP PREDICTOR BY
C COMPARING THE LATTER WITH AN A-B MDEG+2 STEP PREDICTOR.
C IT USES SUBROUTINES PSDIR AND PRDINF(SEE ALSO DOC. ON THESE).
C SELST2 IS WRITTEN BY B. N. LUNDBERG FOR USE WITH THE SOFTWARE ABCON.
C SEE DOCUMENTATION FOR SUBROUTINE PREDIC.
C    1/87
C VERSION 2 : COMPONENTWISE ERROR CONTROL 8/18/87 SEEKS TO KEEP THE
C ERROR IN EACH COMPONENT OF WP LESS THAN THE GIVEN TOLERANCE.
C ** INPUT **
C   MDEG -- DEGREE OF INTERP POLY FOR WHICH STEPSIZE
C           SELECTION FOR PREDICTOR IS SOUGHT(MDEG.LE.MAXDEG-1)
C   N -- DIMENSION OF U IN THE PROBLEM F(U,R) = 0(N.LE.MXDM-1)
C   MAXDEG -- MAXIMUM DEGREE INTERPOLATING POLY ALLOWED.
C   DSP -- STEPSIZE USED IN PREVIOUS STEP.
C   EA  -- ABSOLUTE ERROR REQUIREMENT FOR PREDICTOR
C   ER  -- RELATIVE ERROR REQUIREMENT.
C   HMIN -- MINIMUM STEPSIZE ALLOWED.
C   HMAX -- MAXIMUM STEPSIZE ALLOWED.
C   WC  -- CORRECTED SOLUTION TO F = 0 FROM PREVIOUS STEP.
C   PSYP -- OLD VALUES OF PSY, FOR USE IN PRDINF.
C   PHI  -- SCALED DIVIDED DIFFERENCES FOR USE IN CONSTRUCTING
C           PREDICTOR.
C
C ** OUTPUT **
C    DS -- THE STEPSIZE SELECTED
C    ORFLG-- = .FALSE. , IF ORDER REDUCTION IS NOT INDICATED
C             = .TRUE.  , IF ORDER REDUCTION IS INDICATED.
C
C ** STORAGE AND ARRAYS USED BY OTHER SUBROUTINES CALLED HEREIN **
C     A,B,PSY -- ARRAYS USED IN THE COMPUTATION OF ADAMS-BASHFORTH
C                PREDICTORS(THESE ARRAYS ARE USED IN SUBROUTINE
C                PRDINF AND PSDIR)
C     WP,DIR  -- ARRAYS USED TO STORE THE A-B MDEG+1 STEP PREDICTOR
C                AND PREDICTION DIRECTION COMPUTED IN SUBROUT.PSDIR.
C     G       -- USED IN PRDINF AND PSDIR FOR THE COMPUTATION OF THE
C                A-B PREDICTOR.
C     ALPHA   -- STORES THE STEPSIZE CHANGE FACTOR COMPUTED BELOW.
C     TEMP,WNRM,DIFNRM -- TEMPORARY STORAGE LOCATIONS FOR INTERMEDIATE
C                    RESULTS.
C
C ** SUBROUTINES CALLED **
C    PSDIR  -- GIVES PREDICTED SOLUTION AND PREDICTION DIRECTION.
C    PRDINF -- COMPUTES INFORMATION USED IN PSDIR.
C ** INITIALIZATION **
      DS = DSP
      ORFLG= .FALSE.
C ** COMPUTATION OF NEW DS TO ACHIEVE SPECIFIED ACCURACY **
C     PRINT*,'___________ IN SELST2 ______________________'
         CALL PRDINF(MDEG,MAXDEG,DS,PSYP,PSY,G,A,B)
         CALL PSDIR(MDEG+1,N,MAXDEG,MXDM,DS,WC,WP,DIR,G,B,PHI)
         TEMP =G(MDEG+1,1)*B(MDEG+1)*DS
            WNRM = ABS(WP(1))
            DIFNRM = ABS(TEMP*PHI(1,MDEG+1))
C     PRINT*,'WP1: ',(WP(L),L=1,N+1)
C NUMERATOR AND DENOMINATOR OF ALPHA.
         ANUMER = EA + ER*WNRM
         ADENOM= 2.0*DIFNRM
         IF(ANUMER.GE.ADENOM*(ALIMI**(MDEG+2))) THEN
             ALPHA = ALIMI
           ELSE
              ALPHA = (ANUMER/ADENOM)**(1.0/FLOAT(MDEG+2))
         END IF
         DO 100 L=2,N+1
            WNRM = ABS(WP(L))
            DIFNRM = ABS(TEMP*PHI(L,MDEG+1))
C NUMERATOR AND DENOMINATOR OF ALPHA1.
         ANUMER = EA + ER*WNRM
         ADENOM= 2.0*DIFNRM
         IF(ANUMER.GE.ADENOM*(ALIMI**(MDEG+2))) THEN
             ALPHA1= ALIMI
           ELSE
              ALPHA1= (ANUMER/ADENOM)**(1.0/FLOAT(MDEG+2))
         END IF
         ALPHA = MIN(ALPHA,ALPHA1)
100   CONTINUE
         DS = ALPHA*DS
C        PRINT*,'ALPHA:  ',ALPHA
C FURTHER CHECKS ON DS TO PREPARE FOR LEAVING ROUTINE
510   IF(DS/DSP.LT.ALIMD) THEN
        DS = ALIMD*DSP
        IF(MDEG.GT.0) ORFLG=.TRUE.
      END IF
      IF(DS/DSP.GT.ALIMI) THEN
        DS = ALIMI*DSP
      END IF
      IF(DS.LT.HMIN) THEN
        DS = HMIN
        IF(MDEG.GT.0) ORFLG= .TRUE.
      END IF
550   IF(DS.GT.HMAX) DS = HMAX
                                         RETURN
1000  END
C
C
C
      SUBROUTINE ARCLG2(NP1,MXDM,MAXDEG,MDEG,LEX,DS,TANGNT,DIR,
     +PSYARC,PSY,PSYP,U,PHI,PHIA,B,WC,WCP)
C**************************************************************
C PURPOSE:
C         COMPUTES AN ESIMATE OF THE CHANGE IN ARCLENTH
C         BETWEEN WCP AND WC, FOR USE IN COMPUTING THE
C         NEXT PREDICTION POLYNOMIAL IN AN ADAMS-BASHFORTH
C         METHOD. USES GAUSS-LABOTO QUADRATURE.
C         UPDATES ARRAYS B AND PSY USING THE NEW DS
C ARCLG2 IS WRITTEN BY B. N. LUNDBERG FOR USE WITH THE SOFTWARE ABCON.
C SEE DOCUMENTATION FOR SUBROUTINE PREDIC.
C    SUMMER, 1987
C**************************************************************
      REAL DS,TANGNT(NP1),DIR(NP1),PSYARC(0:MAXDEG),PSY(0:MAXDEG),
     +U(NP1),PHIA(MXDM,0:MAXDEG),
     +B(0:MAXDEG),T,DOT1,DWNRM,TGNRM,PHI(MXDM,0:MAXDEG),PSYP(0:MAXDEG)
     +,WC(NP1),WCP(NP1),DPRIME
      INTEGER NP1,MXDM,MAXDEG,MDEG,LEX,MM
C
C UPDATING MODIFIED DIVIDED DIFFERENCES
C 1: RENORMALIZING : DW/DT*DIR = DIR*DIR-(WC-WCP-2DS*DIR)*DPRIME
C
      MM = MDEG+LEX
C (MM = DEGREE OF INTERPOLATING POLY(BEFORE INTEGRATION) USED TO
C OBTAIN THE PREDICTOR CORRESPONDING TO WC)
      DOT1 = 0.0
      DOT2 = 0.0
      DO 30 L=1,NP1
         DOT1 = DOT1 + TANGNT(L)*DIR(L)
C COMPUTING L TH COMPONENT OF DPRIME AND DOT2 = DW/DT*DIR
         P = PHI(L,MM)*B(MM)
         DO 20 J=MM-1,0,-1
            P = P + PHI(L,J)*B(J)
20       CONTINUE
         DPRIME = (P - DIR(L))/DS
         DOT2 = DOT2+DIR(L)*DIR(L)-(WC(L)-WCP(L)-2.*DS*DIR(L))*DPRIME
30    CONTINUE
      TGNRM = DOT2/DOT1
35    DO 40 L=1,NP1
         PHIA(L,0) = TANGNT(L)*TGNRM
40    CONTINUE
      DO 60 J=1,MM+1
         DO 50 L=1,NP1
            PHIA(L,J) = PHIA(L,J-1)-B(J-1)*PHI(L,J-1)
50       CONTINUE
60    CONTINUE
C ESTIMATING ARCLENGTH INTEGRAL USING GAUSS-LOBATTO RULES.
      IF(MDEG.LE.3) THEN
        SMLBTO = 0.1666666666667D0*(1.0D0 + TGNRM)
        T = 0.5D0
        CALL PDNRM(NP1,MXDM,MAXDEG,MM,DS,T,DWNRM,PSYARC,PSY,PHIA,
     +  U,B)
        SMLBTO = SMLBTO + .66666666666667D0*DWNRM
        GO TO 70
      ELSE IF(MDEG.LE.6) THEN
        SMLBTO = (.83333333333333E-1)*(1.D0+TGNRM)
        T = .27639320225002D0
        CALL PDNRM(NP1,MXDM,MAXDEG,MM,DS,T,DWNRM2,PSYARC,PSY,PHIA,
     +   U,B)
         T = 1.D0 - T
         CALL PDNRM(NP1,MXDM,MAXDEG,MM,DS,T,DWNRM3,PSYARC,PSY,PHIA,
     +   U,B)
         SMLBTO = SMLBTO + (.4166666666667D0)*(DWNRM2+DWNRM3)
         GO TO 70
       ELSE
         SMLBTO =  (.05D0)*(1.D0+TGNRM)
         T = .5D0
         CALL PDNRM(NP1,MXDM,MAXDEG,MM,DS,T,DWNRM3,PSYARC,PSY,PHIA,
     +   U,B)
         SMLBTO = SMLBTO + (.3555555555556D0)*DWNRM3
         T = .172673164646011D0
         CALL PDNRM(NP1,MXDM,MAXDEG,MM,DS,T,DWNRM2,PSYARC,PSY,PHIA,
     +   U,B)
         T = 1.D0- T
         CALL PDNRM(NP1,MXDM,MAXDEG,MM,DS,T,DWNRM4,PSYARC,PSY,PHIA,
     +   U,B)
         SMLBTO = SMLBTO + (.2722222222222D0)*(DWNRM2+DWNRM4)
       END IF
70    DS = DS*SMLBTO
C UPDATING B AND PSY WITH NEW DS
75    B(0) = 1.D0
      PSY(0) = DS
      DO 100 J=1,MDEG+1
         PSY(J) = PSYP(J-1)+DS
         B(J) = B(J-1)*PSY(J-1)/PSYP(J-1)
100   CONTINUE
CPRINTPRINT*,'ESTIMATED ARCLENGTH DS = ',DS
                                                     RETURN
      END
C
      SUBROUTINE PDNRM(NP1,MXDM,MAXDEG,MM,DS,T,DWNRM,PSYARC,
     +PSY,PHIA,U,B)
C******************************************************************
C PURPOSE: USES POLYNOMIAL INTERPOLATION TO ESTIMATE
C THE TANGENT TO THE CURVE AT WC(SN+1 + (T-1)*DS)
C ,SO THAT AS T VARIES FROM 0 TO 1  THE TANGENTS BETWEEN
C WCP AND WC ARE ESTIMATED.
C FOR USE WITH THE ROUTINES ARCLG2 AND PREDIC.
C THE ROUTINE RETURNS THE TANGENT IN U AND ITS NORM IN DWNRM.
C PDNRM IS WRITTEN BY B. N. LUNDBERG FOR USE WITH THE SOFTWARE ABCON.
C SEE DOCUMENTATION FOR SUBROUTINE PREDIC.
C    1/87
C******************************************************************
      REAL DS,T,DWNRM,PSYARC(0:MAXDEG),PSY(0:MAXDEG),PHIA(MXDM,
     +0:MAXDEG),U(NP1),B(0:MAXDEG)
      INTEGER NP1,MXDM,MAXDEG,MM
C
      STEP = (T-1.D0)*DS
      B(0) = 1.0D0
      PSYARC(0) = STEP
      DO 10 J=1,MM+1
         PSYARC(J) = PSY(J-1)+STEP
         B(J) = B(J-1)*PSYARC(J-1)/PSY(J-1)
10    CONTINUE
      DO 20 L=1,NP1
         U(L) = PHIA(L,MM+1)*B(MM+1)
20    CONTINUE
      DO 40 J=MM,0,-1
         DO 30 L=1,NP1
            U(L) = U(L) + B(J)*PHIA(L,J)
30       CONTINUE
40    CONTINUE
      DWNRM = SNRM2(NP1,U,1)
                                              RETURN
      END
      SUBROUTINE COREX(N,MXDM,WC,WP,DIR,TANGNT,D,U,V,STO,STO2,
     +ECOR,ORIENT,DETNEW,KNTROL,NFAC,IPIVOT,ENRM,TOLIT,TOLF,F,DZF,DRF
     +  ,NPRNT)
      PARAMETER(ABERR=0.1,ITMAX=1)
      REAL WC(N+1),WP(N+1),DIR(N+1),TANGNT(N+1),D(MXDM,N+1),
     +     U(N+1),V(N+1),STO(N+1),ECOR,ORIENT,DETNEW,
     +     DIFNRM,ZNRM,TEMP,ENRM,STO2(N+1)
      INTEGER N,MXDM,KNTROL,NFAC,IFLAG,MAXIT,IPIVOT(N+1)
      EXTERNAL F,DZF,DRF
C
C
C
C ** THIS SUBROUTINE CORRECTS A GIVEN GUESS (THE PREDICTED SOLUTION
C    FROM SUBROUTINE PSDIR) BACK TO A SOLUTION TO F = 0, CORRECTING
C    ORTHOGONAL TO THE DIRECTION DIR, USING THE BORDERING ALGORITHM
C    OF H. KELLER.
C    COREX IS WRITTEN BY B. N. LUNDBERG FOR USE WITH THE SOFTWARE ABCON,
C    A FORTRAN SUBROUTINE FOR CONTINUATION ALONG A CURVE DESCRIBED IN:
C    "VARIABLE ORDER ADAMS-BASHFORTH PREDICTORS WITH AN ERROR-
C    STEPSIZE CONTROL FOR CONTINUATION METHODS",  BY B. N. LUNDBERG
C    AND A. B. POORE, SISSC (12), MAY 1991, P.695-723.
C    THE SUBROUTINE ALSO COMPUTES THE DETERMINANT OF THE DERIVATIVE
C    DZF,IN PREPARATION FOR DETECTION OF SINGULAR POINTS
C    ALONG THE CONTINUATION PATH.
C    THE SUBROUTINE ALSO COMPUTES THE NEW UNIT TANGENT
C    DW/DS AT THE NEW POINT (PREDICTED OR CORRECTED).
C    DATE JAN 1987.  UPDATED 5/26/89.
C ** FORTRAN PARAMETERS SPECIFIED ABOVE **
C  ABERR -- AN ABSOLUTE ERROR CONTROL FACTOR USED IN THE STOPPING
C           CRITERION IN THE ITERATION. ABERR =1 GIVES EQUAL RELATIVE
C           AND ABSOLUTE ERROR REQUIREMENTS FOR STOPPING. ABERR.GT.1
C           GIVES MORE WEIGHT TO ABSOLUTE ERROR. ABERR =0.0 GIVES
C           PURE RELATIVE ERROR AND 0.0 .LT. ABERR .LT. 1.0 GIVES
C           SOME ABSOLUTE ERROR CHECK. ABERR = 0.0 IS NOT RECOMMENDED
C           UNLESS ONE KNOWS THE SOLUTION BRANCH IS NEVER NEAR 0.
C  ITMAX -- THE DEFAULT VALUE FOR MAXIMUM NUMBER OF CORRECTOR
C           ITERATIONS.KNTROL MAY BE SET TO ALLOW MORE THAN THIS
C           NUMBER BUT NEVER LESS.
C ** INPUT **
C  MXDM -- ARRAY SIZE PARAMETER;ARRAYS WP,WC,DIR,D,U,V,STO,IPIVOT
C          ARE DIMENSIONED MXDM IN THE CALLING PROGRAM. MXDM IS
C          THE LARGEST VALUE N+1 CAN ASSUME.
C  N -- THE DIMENSION OF U IN THE PROBLEM F(U,R) = 0, THE NUMBER OF EQ'S
C ********************************************************************
C *** CONTROL OPTIONS ***
C  NFAC  =  1 , GIVES DW/DS(TANGENT) AND DET(DZF) AT CORRECTED POINT.
C          ELSE, DW/DS AND DET(DZF) AT PREDICTED VALUE(OR AT
C                PENULTIMATE NEWTON ITERATE IF NEWTONS METHOD IS USED)
C  KNTROL .GT.0 , IF CHORD METHOD(ONE FACTORIZATION) IS REQUESTED.
C         .LT.0 , IF NEWTONS METHOD IS REQUESTED.
C         .GT.ITMAX OR .LT. -ITMAX , IF MAXIT = ABS(KNTROL).GT.
C                            ITMAX ITERATIONS ARE TO BE ALLOWED.
C ********************************************************************
C  WP -- PREDICTOR OF THE SOLUTION AS ON OUTPUT FROM SUBROUTINE PSDIR.
C  DIR -- PREDICTION DIRECTION AS ON OUTPUT FROM SUBROUTINE PSDIR.
C  ECOR --  CORRECTOR ERROR SPECIFICATION. USED TO DETERMINE THE
C           STOPPING CRITERION FOR THE CHORD OR NEWTON ITERATIONS.
C  ORIENT --  = 1.0  IF POSITIVE ORIENTATION(CHANGE IN PARAMETER R)
C         = -1.0 IF NEGATIVE ORIENTATION IS DESIRED.
C
C ** OUTPUT **
C  WC -- THE CORRECTED SOLUTION
C  TANGNT  -- THE NEW RIGHT HAND SIDE OF THE DAVIDENKO ODE
C           WITH ARCLENGTH PARAMETERIZATION.  SEE DOC. ON NFAC ABOVE
C           FOR OPTIONS.
C  DETNEW -- THE DETERMINANT OF DZF, THE DERIVATIVE OF F WITH RESPECT
C            TO THE FIRST N VARIABLES(AT THE PREDICTED OR CORRECTED
C            SOLUTION POINT AS SPECIFIED BY NFAC).
C  KNTROL  -- RETURNS THE NUMBER OF CORRECTOR STEPS TAKEN.
C             OR:*** IS SET TO 0(ZERO) IF CORRECTOR FAILS ***
C
C  ENRM  -- THE EUCLIDEAN NORM OF W1-WP = W1-W0
C ** STORAGE **
C  D(N,N) (FOR DERIVATIVE),U,V,STO,IPIVOT
C  DIFNRM,ZNRM -- USED IN STORING RESULTS RELATED TO STOPPING
C                 CRITERION.
C  TEMP -- A TEMPORARY STORAGE LOCATION FOR VARIOUS INTERMEDIATE
C          RESULTS.
C  IFLAG -- USED IN CALLS TO SGEFA AND AS AN INDICATOR VARIABLE
C           IN THE STOPPING CRITERION.
C  DELR -- A STORAGE VARIABLE USED IN COMPUTING THE CORRECTOR ITERATES.
C  MAXIT -- STORES THE LIMIT ON THE NUMBER OF CORRECTOR ITERATIONS.
C            (MAXIT IS AT LEAST ITMAX)
C
      MAXIT=MAX(ABS(KNTROL),ITMAX)
      NP1 = N + 1
      DO 10 L=1,NP1
         WC(L) = WP(L)
 10   CONTINUE
      CALL DZF(N,MXDM,D,WC)
      CALL SGEFA(D,MXDM,N,IPIVOT,IFLAG)
      IF(IFLAG.GT.0) THEN
        IF(NPRNT.GE.3) THEN
          PRINT*,'ZERO PIVOT IN COREX  PLU'
        END IF
        IRCODE = -4
        KNTROL = 0
                                       RETURN
      END IF
      IF(NFAC.NE.1) THEN
        CALL SGEDI(D,MXDM,N,IPIVOT,DETNEW,V,10)
      END IF
         CALL DRF(N,U,WC)
         CALL SGESL(D,MXDM,N,IPIVOT,U,0)
         CALL F(N,V,WC)
         FNRMP = ABS(V(N))
         DO 12 L=1,N-1
            FNRMP = MAX(FNRMP,ABS(V(L)))
12       CONTINUE
        F0 = FNRMP
CPRINT  PRINT*,'**** IN COREX ****'
CPRINT  PRINT*,'FNRM(WP) = ',FNRMP
19    DO 50 I=1,MAXIT
         CALL SGESL(D,MXDM,N,IPIVOT,V,0)
         TEMP = DIR(NP1)
         DELR = (WC(NP1)-WP(NP1))*TEMP
         DO 20 L=1,N
            DELR = DELR + (WC(L)-WP(L)-V(L))*DIR(L)
            TEMP = TEMP - U(L)*DIR(L)
  20     CONTINUE
         IFLAG =0
         DELR = -DELR/TEMP
         DIFNRM = ABS(DELR)
         DWNRM = DIFNRM
         WC(NP1)= WC(NP1)+ DELR
         ZNRM =(ABS(WC(NP1))+ ABERR)*ECOR
         IF(DIFNRM.GT.ZNRM) IFLAG = 1
         DO 30 L=1,N
            TEMP= -V(L) - DELR*U(L)
            WC(L) = WC(L) + TEMP
            DIFNRM= ABS(TEMP)
            DWNRM = MAX(DWNRM,DIFNRM)
            ZNRM = (ABS(WC(L))+ABERR)*ECOR
            IF(DIFNRM.GT.ZNRM) IFLAG = 1
  30     CONTINUE
C ** TESTING FOR CONVERGENCE AND NON-CONVERGENCE **
         CALL F(N,V,WC)
         FNRM = ABS(V(1))
         DO 35 L=2,N
            FNRM = MAX(ABS(V(L)),FNRM)
35       CONTINUE
         IF((IFLAG.LE.0).AND.(FNRM.LT.ECOR*MAX(F0,1.))) GO TO 60
CPRINT   PRINT*,'FNRM = ',FNRM,'  DWNRM = ',DWNRM
         IF(FNRM.GT.TOLF*FNRMP) THEN
           KNTROL = 0
           IF(NPRNT.GT.3) THEN
           PRINT*,'FNRM/FNRMP TO LARGE IN COREX'
           END IF
           GO TO 1000
         END IF
         IF(I.GT.1) THEN
           IF(DWNRM.GT.TOLIT*DWNRMP) THEN
             KNTROL = 0
             IF(NPRNT.GT.3) THEN
             PRINT*,'DWNRM/DWNRMP TO LARGE'
             END IF
             GO TO 1000
           END IF
         END IF
         FNRMP = FNRM
         DWNRMP = DWNRM
C ** NEWTONS METHOD OPTION **
         IF(KNTROL.LT.0) THEN
           CALL DZF(N,MXDM,D,WC)
           CALL SGEFA(D,MXDM,N,IPIVOT,IFLAG)
           IF(IFLAG.GT.0) THEN
             IF(NPRNT.GE.3) THEN
               PRINT*,'ZERO PIVOT IN COREX  PLU'
             END IF
             IRCODE = -4
             KNTROL = 0
                                       RETURN
           END IF
           CALL DRF(N,U,WC)
           CALL SGESL(D,MXDM,N,IPIVOT,U,0)
         END IF
 50   CONTINUE
C * INDICATING CONVERGENCE NOT OBTAINED AND EXITING ROUTINE *
      KNTROL = 0
           IF(NPRNT.GT.3) THEN
           PRINT*,'ITERATION LIMIT EXCEEDED IN CORRECTOR'
           END IF
                                          RETURN
C INDICATING NUMBER OF ITERATIONS IN CASE CONVERGENCE OBTAINED
60    KNTROL = I
C COMPUTING ERROR
      ENRM = 0.
      DO 62 L=1,NP1
         TEMP = ABS(WC(L) -WP(L))
         ENRM = MAX(ENRM,TEMP)
62    CONTINUE
C OPTION NFAC = 1: TANGENT AND DETERMINANT AT CORRECTED VALUE.
      IF(NFAC.EQ.1) THEN
        CALL DZF(N,MXDM,D,WC)
        CALL SGEFA(D,MXDM,N,IPIVOT,IFLAG)
           IF(IFLAG.GT.0) THEN
             IF(NPRNT.GE.3) THEN
               PRINT*,'ZERO PIVOT IN COREX  PLU'
             END IF
             IRCODE = -4
             KNTROL = 0
                                       RETURN
         END IF
        CALL DRF(N,U,WC)
C COMPUTING NEW DETERMINANT DET(DZF)
        CALL SGEDI(D,MXDM,N,IPIVOT,DETNEW,V,10)
C COMPUTING NEW DW/DS TO SET UP FOR NEXT PREDICTOR STEP.
      CALL DWDS(D,N,MXDM,U,TANGNT,IPIVOT,ORIENT)
      ELSE
      TEMP = 1.0
      DO 150 L=1,N
         TEMP = TEMP + U(L)**2
150    CONTINUE
      TEMP = ORIENT/SQRT(TEMP)
      DO 200 L=1,N
         TANGNT(L) = -U(L)*TEMP
200   CONTINUE
      TANGNT(N+1) = TEMP
      END IF
C
C
C
1000                                                RETURN
      END
C
C
C
C *************  D W D S ********************************************
C
C
      SUBROUTINE DWDS(DZ,N,MXDM,DR,TANGNT,IPIVOT,ORIENT)
      REAL DZ(MXDM,N+1),DR(N+1),TANGNT(N+1),ORIENT,TEMP
      INTEGER N,MXDM,IPIVOT(N+1)
C THIS ROUTINE COMPUTES DW/DS, THE RIGHT HAND SIDE OF THE DAVIDENKO
C O.D.E. WITH ARCLENGTH PARAMETERIZATION. DW/DS IS THE UNIT TANGENT
C TO THE SOLUTION BRANCH THE CONTINUATION CURVE AND IS USED TO COMPUTE
C THE PREDICTOR IN PREDICTOR-CORRECTOR , ARCLENGTH CONTINUATION.
C DWDS IS WRITTEN BY B. N. LUNDBERG FOR USE WITH THE SOFTWARE ABCON.
C SEE DOCUMENTATION FOR SUBROUTINE PREDIC.
C   UPDATED 5/26/89
C
C ** INPUT **
C  DZ(N,N) -- CONTAINS THE LU FACTORIZATION OF THE MATRIX OF
C             DF(U,R)/DU AS ON OUTPUT FROM SUBROUTINES DZF
C             AND SGEFA.
C  DR(N)   -- CONTAINS THE DERIVATIVE DF(U,R)/DR AS ON OUTPUT
C             FROM SUBROUTINE DRF.
C  N --       THE SIZE OF THE VECTOR U IN THE PROBLEM F(U,R) = 0.
C  MXDM  --   THE MAX DIMENSION OF (U,R) IN THE PROBLEM F(U,R)=0.
C             D IS DIMENSIONED MXDM BY MXDM IN THE CALLING PROG.
C
C  IPIVOT --  THE VECTOR WHICH, ON OUTPUT FROM SGEFA, RECORDS
C             THE PIVOTS FOR USE IN A CALL OF SUBROUTINE SGESL.
C  ORIENT --  THE ORIENTATION OF THE CONTINUATION.
C             = 1.0 FOR POSITIVE MOVEMENT IN R
C              -1.0 FOR NEGATIVE MOVEMENT IN R.
C ** OUTPUT **
C  TANGNT(N+1) -- THE VECTOR DW/DS.
C ** STORAGE **
C  TEMP -- USED FOR STORING AND COMPUTING A VECTOR NORM.
C
      CALL SGESL(DZ,MXDM,N,IPIVOT,DR,0)
      TEMP = 1.0
      DO 50 L=1,N
         TEMP = TEMP + DR(L)**2
50    CONTINUE
      TEMP = ORIENT/SQRT(TEMP)
      DO 100 L=1,N
         TANGNT(L) = -DR(L)*TEMP
100   CONTINUE
      TANGNT(N+1) = TEMP
                                        RETURN
      END
C
      SUBROUTINE SINGCK(N,MXDM,WC,WCP,D,TANGNT,TNGNTP,V,DZF,DRF,INDIC,
     + IPIVOT)
C  PURPOSE: THIS SUBROUTINE CHECKS THE DETERMINANTS OF THE 
C LARGER SYSTEM MATRICES AT CURRENT AND PREVIOUS POINTS TO DETERMINE
C WHETHER A SIMPLE FOLD OR A SIMPLE BIFURCATION POINT HAS BEEN PASSED.
C THIS ROUTINE USES PLU ROUTINE SGEFA TO ALLOW
C EASY COMPUTATION OF DETERMINANTS OF AN (N+1)x(N+1) MATRIX
C SINGCK IS WRITTEN BY B. N. LUNDBERG FOR USE WITH THE SOFTWARE ABCON.
C SEE DOCUMENTATION FOR SUBROUTINE PREDIC.
C    REWRITE: 6/90, 5/93.
C   INDIC = 1, IF A FOLD POINT IS INDICATED.
C         = 2, IF A BIFURCATION POINT IS INDICATED.
C
      REAL WC(N+1), WCP(N+1),TANGNT(N+1),TNGNTP(N+1),DETS1,DETS2,
     &DETL1,DETL2,D(MXDM,N+1),V(N+1)
      INTEGER INDIC, IPIVOT(N+1),IFLAG
      EXTERNAL DZF,DRF
C
C
      CALL DZF(N,MXDM,D,WCP)
      CALL DRF(N,V,WCP)
      CALL MAKMAT(N,MXDM,D,V,TNGNTP)
C
CPRINT        PRINT*,'**********  PREVIOUS JACOBIAN *********'
C             DO 278 LL=1,N+1
C                PRINT*,(D(LL,II),II=1,N+1)
C278           CONTINUE
C           PRINT*,'*********************************'
      IFLAG = 0
      CALL SGEFA(D,MXDM,N+1,IPIVOT,IFLAG)
      CALL SGEDI(D,MXDM,N+1,IPIVOT,DETL1,V,10)
C
C     PRINT*,'PREVIOUS DETS: SMALL = ',DETS1,' LARGE = ',DETL1
C COMPUTING NEW DETERMINANTS
C
C
      CALL DZF(N,MXDM,D,WC)
      CALL DRF(N,V,WC)
      CALL MAKMAT(N,MXDM,D,V,TANGNT)
C
C           PRINT*,'********** CURRENT JACOBIAN *********'
C             DO 288 LL=1,N+1
C                PRINT*,(D(LL,II),II=1,N+1)
C288           CONTINUE
C           PRINT*,'*********************************'
      CALL SGEFA(D,MXDM,N+1,IPIVOT,IFLAG)
      CALL SGEDI(D,MXDM,N+1,IPIVOT,DETL2,V,10)
C     PRINT*,'CURRENT DETS: SMALL = ',DETS2,' LARGE = ',DETL2
C TESTING FOR FOLD OR BIFURCATION
C
      INDIC = 1
      IF(DETL1*DETL2.LE.0.) THEN
        INDIC = 2
C       PRINT*,'BIFURCATION POINT DETECTED IN SINGCHK'
      ELSE
C       PRINT*,'FOLD POINT DETECTED IN SINGCHK'
      END IF
                            RETURN
      END
C
C
C
      SUBROUTINE MAKMAT(N,MXDM,D,DR,TAN)
      REAL D(MXDM,N+1),DR(N),TAN(N+1)
      INTEGER N,MXDM
C CONSTRUCTS LARGER SYSTEM MATRIX FOR ROUTINE SINGCK
C MAKMAT IS WRITTEN BY B. N. LUNDBERG FOR USE WITH THE SOFTWARE ABCON.
C SEE DOCUMENTATION FOR SUBROUTINE PREDIC OR ABCON.
C    REWRITE 6/90
C  
      DO 20 I=1,N
         D(N+1,I) = TAN(I)
         D(I,N+1) = DR(I)
20    CONTINUE
      D(N+1,N+1) = TAN(N+1)
                                RETURN
      END
C
C *************  F I N A L P ****************************************
C
C
      SUBROUTINE FINALP(N,MDEG,MXDM,MAXDEG,D,WCP,WP,STO1,STO2,DS,
     +DIR,EFIN,RLOW,RHIGH,LEX,KNTROL,G,PSYP,PSY,A,B,PHI,IPIVOT,WC,
     &NCOMP,VCOMP,F,DZF,DRF,NPRNT,FINAL,IFIN,EC,IRCODE)
      PARAMETER(MXIT=100,TOLSTP=1.E-8,ABERR=0.1)
      REAL D(MXDM,N+1),WCP(N+1),WP(N+1),STO1(N+1),
     +STO2(N+1),DS,DIR(N+1),EFIN,RLOW,RHIGH,PSYP(0:MAXDEG),
     +PSY(0:MAXDEG),G(0:MAXDEG,1:MAXDEG),A(0:MAXDEG),EC,
     +B(0:MAXDEG),PHI(MXDM,0:MAXDEG),ZNRM,DIFNRM,WC(N+1)
      INTEGER N,MDEG,MXDM,MAXDEG,LEX,KNTROL,IPIVOT(N+1),IFLAG
     +        ,IFIN,IRCODE
      EXTERNAL F,DZF,DRF
      LOGICAL FINAL
C
C
C *** PURPOSE ***
C     THIS SUBROUTINE CHECKS WHETHER OR NOT GIVEN PARAMETER 
C     BOUNDS OR TARGET VALUE HAVE BEEN CROSSED.  IF YES,
C     PREPARATION IS MADE FOR THE FINAL STEP OF THE CONTINUATION.
C     A PREDICTION IS MADE WHOSE ICOMP COMPONENT MATCHES THE
C     TARGET VALUE VCOMP.  ICOMP IS + OR - NCOMP, OR N+1,
C     VCOMP IS VTARG, RHIGH OR RLOW.
C     THIS SUBROUTINE THEN COMPUTES A SOLUTION TO F(U,R) = 0
C     WHERE THE ICOMP COMPONENT OF W=(U,R) IS HELD FIXED.
C     IT USED THE PRESENT ORDER A-B PREDICTOR TO PREDICT
C     A POINT ON THE CURVE,  AND THEN CORRECTS
C     THIS GUESS USING NEWTON(EFIN.LT.0) OR CHORD(EFIN.GT.0)
C     METHOD.   THE STOPPING ACCURACY REQUIREMENT IS ABS(EFIN).
C     THE SOLUTION IS RETURNED IN WC,  THE PREDICTOR IN WP.
C     NO CORRECTION IS PERFORMED IN THIS ROUTINE IF EFIN = 0.
C     OF IFIN.LE.0, THE ROUTINE PREPARES FOR A CORRECTION
C     USING THE USUAL CORRECTOR ROUTINE COREX, WITH TARGET
C     PARAMETER FIXED (IFIN.LT.0) OR ITERATES RESTRICTED TO
C     HYPERPLANE ORTHOGONAL TO PREDICTOR DIRECTION.
C
C  FINALP IS WRITTEN BY B. N. LUNDBERG FOR USE WITH THE SOFTWARE ABCON,
C  A FORTRAN SUBROUTINE FOR CONTINUATION ALONG A CURVE DESCRIBED IN:
C  "VARIABLE ORDER ADAMS-BASHFORTH PREDICTORS WITH AN ERROR-
C  STEPSIZE CONTROL FOR CONTINUATION METHODS",  BY B. N. LUNDBERG
C  AND A. B. POORE, SISSC (12), MAY 1991, P.695-723.
C  5/26/89:   UPDATE: 7/18/92. 
      NP1 = N + 1
C ****** TESTING FOR SPECIAL POINT TERMINATION  *************
C
C CHECKING WHETHER PREDICTION OR PREVIOUS POINT HAS CROSSED BOUNDARY
C
      IF(NCOMP.GT.0) THEN
        IF((WCP(NCOMP).GE.VCOMP).OR.(WP(NCOMP).GE.VCOMP)) THEN
          IF(NPRNT.GT.3) THEN
          PRINT*,'UPPER PARAMETER LIMIT PASSED ',VCOMP
          END IF
          FINAL = .TRUE.
          ICOMP = NCOMP
          VTARG = VCOMP
        END IF
      ELSE IF(NCOMP.LT.0) THEN
        IF((WCP(-NCOMP).LE.VCOMP).OR.(WP(-NCOMP).LE.VCOMP)) THEN
          IF(NPRNT.GT.3) THEN
             PRINT*,'LOWER PARAMETER LIMIT PASSED ',VCOMP
          END IF
          FINAL = .TRUE.
          ICOMP = -NCOMP
          VTARG = VCOMP
        END IF
      ELSE IF(NCOMP.EQ.0) THEN
        IF((WCP(NP1).GE.RHIGH).OR.(WP(NP1).GE.RHIGH)) THEN
          FINAL = .TRUE.
          ICOMP = NP1
          VTARG = RHIGH
        END IF
        IF((WCP(-NCOMP).LE.VCOMP).OR.(WP(-NCOMP).LE.VCOMP)) THEN
          FINAL = .TRUE.
          ICOMP = NP1
          VTARG = RLOW
        END IF
      END IF
C
C *** RETURN IF ALL CHECKS ARE FALSE, RETURN, CONTINUE TRACKING ***
      IF(.NOT.FINAL)                          RETURN
C
C *** P R E D I C T O R  C O M P U T A T I O N ********************
C
C ITERATION TO DETERMINE STEPSIZE TO REACH (U,R)
      DIR(ICOMP) = (WP(ICOMP) - WCP(ICOMP))/DS
      DO 200 I=1,MXIT
         DS = (VTARG - WCP(ICOMP))/DIR(ICOMP)
         CALL PRDINF(MDEG,MAXDEG,DS,PSYP,PSY,G,A,B)
         CALL PSDIR(MDEG+LEX,N,MAXDEG,MXDM,DS,WCP,WP,DIR,G,B,PHI)
         IF(ABS(VTARG-WP(ICOMP)).LT.TOLSTP) GO TO 210
200   CONTINUE
      IF(NPRNT.GT.3) THEN
      PRINT*,'FAILURE IN STEP SELECTION FOR FINAL ITERATION'
      END IF
      IRCODE = -2
      KNTROL = 0
                                          RETURN
210   CONTINUE
C
C ***  PREPARING FOR CORRECTION OR RETURN IN MAIN  ****
C
C IF EFIN = 0, NO CORRECTION TO FINAL POINT
      IF(EFIN.EQ.0.) THEN
        DO 220 L=1,NP1
          WC(L) = WP(L)
220     CONTINUE
                   RETURN
      END IF
C
C IF IFIN < 0, PREPARE FOR FINAL CORRECTION IN COREX
C
      IF(IFIN.LT.0.) THEN
C ADJUSTING PREDICTOR DIRECTION FOR CORRECTION WITH TARGET
C PARAMETER FIXED ('VERTICAL CORRECTION')
        DO 230 L=1,N+1
           DIR(L) = 0.
230     CONTINUE
        DIR(ICOMP) = 1.
        WP(ICOMP) = VTARG
        EC = ABS(EFIN)
        KNTROL = MAX(ABS(KNTROL),ABS(IFIN))*SIGN(1.E0,EFIN)
      END IF
      IF(IFIN.LE.0.)                 RETURN
C
C ****** F I N A L  C O R R E C T I O N ***********************
C
C IF IFIN > 0, NEWTON'S METHOD CORRECTOR (NO BORDERING ALGORITHM)
C WITH ICOMP COMPONENT KEPT FIXED AT VTARG
      DO 235 L=1,N+1
         WC(L) = WP(L)
235   CONTINUE
      WC(ICOMP) = VTARG
      ITMAX = MAX(ABS(KNTROL),IFIN)
      DO 300 I=1,ITMAX
         JFLAG = 0
         IF((EFIN.LT.0.).OR.(I.EQ.1)) THEN
C NEWTONS METHOD OPTION AND FIRST STEP
           CALL DZF(N,MXDM,D,WC)
           IF(ICOMP.LT.NP1) THEN
             DO 240 J=ICOMP,N-1
                DO 250 L=1,N
                   D(L,J) = D(L,J+1)
250             CONTINUE
240          CONTINUE
             CALL DRF(N,STO1,WC)
             DO 270 L=1,N
                D(L,N) = STO1(L)
270          CONTINUE
           END IF
           IFLAG = 1
           CALL SGEFA(D,MXDM,N,IPIVOT,IFLAG)
           IF(IFLAG.GT.0) THEN
             IF(NPRNT.GT.3) THEN
               PRINT*,'ZERO PIVOT IN PLU IN FINAL********'
             END IF
             IRCODE = -4
             KNTROL = 0
             DO 275 L=1,NP1
               WC(L) = WCP(L)
275          CONTINUE
                                                RETURN
           END IF
         END IF
         CALL F(N,STO2,WC)
         CALL SGESL(D,MXDM,N,IPIVOT,STO2,0)
C UPDATE ALL BUT ICOMP ENTRY (WC(ICOMP) FIXED AT VTARG)
         DO 280 L=1,ICOMP-1
            WC(L) = WC(L) - STO2(L)
            DIFNRM = ABS(STO2(L))
            ZNRM = (ABS(WC(L)) + ABERR)*ABS(EFIN)
            IF(DIFNRM.GT.ZNRM) JFLAG = 1
280      CONTINUE
         DO 290 L=ICOMP,N
            WC(L+1) = WC(L+1) - STO2(L)
            DIFNRM = ABS(STO2(L))
            ZNRM = (ABS(WC(L+1)) + ABERR)*ABS(EFIN)
            IF(DIFNRM.GT.ZNRM) JFLAG = 1
290      CONTINUE
C TESTING FOR CONVERGENCE
         IF(JFLAG.EQ.0)  GO TO 310
300   CONTINUE
C INDICATING FAILURE TO CONVERGE
      KNTROL = 0
      IRCODE = -3
      DO 305 L=1,NP1
         WC(L) = WCP(L)
305   CONTINUE
                                       RETURN
310   CONTINUE
C INDICATING NUMBER OF ITERATIONS IN FINAL CORRECTION
      KNTROL = I
                                       RETURN
      END
C
C 
C 
C 
C*********** FUNCTIONS AND DERIVATIVES ******************* 
C 
      SUBROUTINE F1(N,RESULT,ARG)
      REAL RESULT(N),ARG(N+1) 
      INTEGER N 
      COMMON KJAC,KFUNC 
C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION 
C     WE WISH TO CONTINUE ALONG. ARG CONTAINS THE VALUE OF (X(N),R) 
C     AT WHICH F IS COMPUTED.  RESULT(N)  CONTAINS F(X,R) ON RETURN.
C     ARG,N ARE UNCHANGED.
C 
C 
C ** SCHWETLICK EXAMPLE #3 ************ 
C 
C 
      KFUNC = KFUNC + 1 
      RESULT(1) = 200.*ARG(1)**3-200.*ARG(1)*ARG(2)+ARG(1)-1. 
      RESULT(2) = -100.*ARG(1)**2+110.1*ARG(2)+9.9*ARG(4)-20. 
      RESULT(3) = 180.*ARG(3)**3-180.*ARG(3)*ARG(4)+ARG(3)-1. 
      RESULT(4) = 9.9*ARG(2)-90.*ARG(3)**2+100.1*ARG(4)-20. 
C CONSTRUCTING STANDARD HOMOTOPY
      RESULT(1) = RESULT(1) - (1.-ARG(5))*(-6004.)
      RESULT(2) = RESULT(2) - (1.-ARG(5))*(-1040.)
      RESULT(3) = RESULT(3) - (1.-ARG(5))*(-5404.)
      RESULT(4) = RESULT(4) - (1.-ARG(5))*(-940.) 
                                      RETURN
      END 
C 
C 
C 
      SUBROUTINE DZF1(N,MXDM,DZ,ARG) 
      REAL DZ(MXDM,N), ARG(N+1) 
      INTEGER N,MXDM
      COMMON KJAC,KFUNC 
C ***  COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F
C      WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED
C      IN ARG(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N).
C 
C ********* SCHWETLICK EXAMPLE # 3 ************** 
C 
      KJAC = KJAC+1 
C 
      DO 20 J=1,N 
         DO 10 I=1,N
            DZ(I,J) = 0.
10       CONTINUE 
20    CONTINUE
C 
      DZ(1,1) =600.*ARG(1)**2-200.*ARG(2)+1.
      DZ(1,2) = -200.*ARG(1)
C 2ND COMPONENT 
      DZ(2,1) = -200.*ARG(1)
      DZ(2,2) = 110.1 
      DZ(2,4) = 9.9 
C THIRD COMPONENT 
      DZ(3,3) = 540.*ARG(3)**2-180.*ARG(4)+1. 
      DZ(3,4) = -180.*ARG(3)
C FOURTH COMPONENT
      DZ(4,2) = 9.9 
      DZ(4,3) = -180.*ARG(3)
      DZ(4,4) = 100.1 
                                     RETURN 
      END 
C 
C 
C 
      SUBROUTINE DRF1(N,DR,ARG)
      REAL DR(N),ARG(N+1) 
      INTEGER N 
C ***  COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F
C      WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN ARG(N+1)
C      THE DERIVATIVE IS RETURNED IN DR(N). 
C 
      DR(1) = -6004.
      DR(2) = -1040.
      DR(3) = -5404.
      DR(4) = - 940.
                                  RETURN
      END 
C FUNCTION SCHWETLICK # 4 ******************************************* 
C 
C **** FUNCTION  AND DERIVATIVES **** 
C 
C 
      SUBROUTINE F2(N,RESULT,ARG)
      REAL RESULT(N),ARG(N+1) 
      INTEGER N 
      COMMON KJAC,KFUNC 
C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION 
C     WE WISH TO CONTINUE ALONG. ARG CONTAINS THE VALUE OF (X(N),R) 
C     AT WHICH F IS COMPUTED.  RESULT(N)  CONTAINS F(X,R) ON RETURN.
C     ARG,N ARE UNCHANGED.
C 
C 
C ** SCHWETLICK EXAMPLE #4 ************ 
C 
      FI(Y) = 5.6E-8*(EXP(25.*Y) - 1.)
      FU(Y) = 7.65*ATAN(1962.*Y)
C 
      KFUNC = KFUNC + 1 
      RESULT(1) = 1.E-4*(ARG(1) - ARG(3))+(ARG(1)-ARG(2))/39.+
     1            (ARG(1)+ARG(7))/51. 
      RESULT(2) = .1*(ARG(2)-ARG(6))+(ARG(2)-ARG(1))/39.+FI(ARG(2)) 
      RESULT(3) = 1.E-4*(ARG(3)-ARG(1))+(ARG(3)-ARG(4))/25.5
      RESULT(4) = (ARG(4)-ARG(3))/25.5+ARG(4)/0.62+ARG(4)-ARG(5)
      RESULT(5) = (ARG(5)-ARG(6))/13.+ARG(5)-ARG(4)+FI(ARG(5))
      RESULT(6) = (ARG(6)-ARG(2))*.1+(ARG(6)-ARG(5))/13.+ARG(6) 
     1           -FU(ARG(3)-ARG(1))/0.201 
                                      RETURN
      END 
C 
C 
C 
      SUBROUTINE DZF2(N,MXDM,DZ,ARG) 
      REAL DZ(MXDM,N), ARG(N+1),FN(1001)
      INTEGER N,MXDM
      COMMON KJAC,KFUNC 
C ***  COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F
C      WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED
C      IN ARG(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N).
C 
C ********* SCHWETLICK EXAMPLE # 4 ************** 
C 
      DI(Y) = 25.*5.6E-8*EXP(25.*Y) 
      DU(Y) = 1962.*7.65/(1.+(1962.*Y)**2)
      KJAC = KJAC+1 
C 
      DO 20 J=1,N 
         DO 10 I=1,N
            DZ(I,J) = 0.
10       CONTINUE 
20    CONTINUE
C 
      DZ(1,1) =1.E-4+1./39.+1./51.
      DZ(1,2) = -1./39. 
      DZ(1,3) = -1.E-4
C 2ND COMPONENT 
      DZ(2,1) = -1./39. 
      DZ(2,2) = .1+1./39. + DI(ARG(2))
      DZ(2,6) = -.1 
C THIRD COMPONENT 
      DZ(3,1) = -1.E-4
      DZ(3,3) = 1.E-4+1./25.5 
      DZ(3,4) = -1./25.5
C FOURTH COMPONENT
      DZ(4,3) = -1./25.5
      DZ(4,4) = 1./25.5+1./0.62+1.
      DZ(4,5) = -1. 
C FIFTH COMPONENT 
      DZ(5,4) = -1. 
      DZ(5,5) = 1./13.+1.+DI(ARG(5))
      DZ(5,6) = -1./13. 
C SIXTH COMPONENT 
      DZ(6,1) = DU(ARG(3) - ARG(1))/0.201 
      DZ(6,2) = -.1 
      DZ(6,3) = -DU(ARG(3)-ARG(1))/0.201
      DZ(6,5) = -1./13. 
      DZ(6,6) = .1+1./13.+1.
                                     RETURN 
      END 
C 
C 
C 
      SUBROUTINE DRF2(N,DR,ARG)
      REAL DR(N),ARG(N+1) 
      INTEGER N 
C ***  COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F
C      WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN ARG(N+1)
C      THE DERIVATIVE IS RETURNED IN DR(N). 
C 
      DO 10 I=1,N 
         DR(I) = 0.0
10    CONTINUE
      DR(1) = 1./51.
                                       RETURN 
      END 
C 
C 
C************ F & R FUNCTION WITH STANDARD HOMOTOPY 
C 
C **** FUNCTION  AND DERIVATIVES **** 
C 
      SUBROUTINE F3(N,RESULT,ARG)
C ** FREUDENSTIEN AND ROTH ** 
      REAL RESULT(N),ARG(N+1) 
      INTEGER N 
      COMMON KJAC,KFUNC 
      KFUNC = KFUNC+1 
C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION 
C     WE WISH TO CONTINUE ALONG. ARG CONTAINS THE VALUE OF (X(N),R) 
C     AT WHICH F IS COMPUTED.  RESULT(N)  CONTAINS F(X,R) ON RETURN.
C     ARG,N ARE UNCHANGED.
C 
C 
      RESULT(1)=ARG(1)+5.*ARG(2)**2-ARG(2)**3-2.*ARG(2) 
     1  -47.+ 34.*ARG(3)
      RESULT(2)=ARG(1)+ARG(2)**2+ARG(2)**3-14.*ARG(2) 
     1   -39. +10.*ARG(3) 
                                        RETURN
      END 
C 
C 
      SUBROUTINE DZF3(N,MXDM,DZ,ARG) 
      REAL DZ(MXDM,N), ARG(N+1) 
      INTEGER N,MXDM
      COMMON KJAC,KFUNC 
      KJAC = KJAC + 1 
C ***  COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F
C      WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED
C      IN ARG(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N).
C 
      DZ(1,1) =  1.0
      DZ(2,1) = 1.0 
      DZ(1,2) = -3.*ARG(2)**2+10.*ARG(2)-2. 
      DZ(2,2) = 3.*ARG(2)**2+2.*ARG(2)-14.
                                       RETURN 
      END 
C 
C 
C 
      SUBROUTINE DRF3(N,DR,ARG)
      REAL DR(N),ARG(N+1) 
      INTEGER N 
C ***  COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F
C      WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN ARG(N+1)
C      THE DERIVATIVE IS RETURNED IN DR(N). 
C 
      DR(1) = 34. 
      DR(2) = 10. 
                                       RETURN 
      END 
C ***** FUNCTION 1 FROM WATSON *****************
C 
C 
      SUBROUTINE F4(NV,RESULT,ARG) 
      REAL RESULT(NV),ARG(NV+1) 
      INTEGER NV
      COMMON KJAC,KFUNC 
C ** WATSON 1 **
C 
C 
      KFUNC = KFUNC + 1 
      SUM = 0.
      DO 10 L=1,NV
         SUM = SUM + ARG(L)**3
10    CONTINUE
      TEMP = 1./(2.*FLOAT(NV))
      DO 20 L=1,NV
         RESULT(L) = (SUM + FLOAT(L))*TEMP
20    CONTINUE
C FORMING HOMOTOPY (REG.) 
      TEMP = ARG(NV+1)
      DO 30 L=1,NV
         RESULT(L) = TEMP*(ARG(L)-RESULT(L))+ (1.-TEMP)*ARG(L)
30    CONTINUE
                                        RETURN
      END 
C 
C 
      SUBROUTINE DZF4(NV,MXDM,DZ,ARG)
      REAL DZ(MXDM,NV), ARG(NV+1),FN(101) 
      INTEGER N,NV
      COMMON KJAC,KFUNC 
      KJAC = KJAC+1 
      SUM = 0.
       N = NV 
       TEMP = 1./(2.*FLOAT(N))
      TEMP =-3.*TEMP*ARG(N+1) 
      DO 40 I=1,N 
         TEMP1 = TEMP*ARG(I)**2 
         DO 30 L=1,N
            DZ(L,I) = TEMP1 
30       CONTINUE 
         DZ(I,I) = 1.  + DZ(I,I)
40     CONTINUE 
                        RETURN
      END 
      SUBROUTINE DRF4(N,DR,ARG)
      REAL DR(N),ARG(N+1) 
      INTEGER N 
C ***  COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F
C      WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN ARG(N+1)
C      THE DERIVATIVE IS RETURNED IN DR(N). 
C 
      TEMP = 1./(2.*FLOAT(N)) 
      SUM = 0.
      DO 10 L=1,N 
         SUM = SUM + ARG(L)**3
10    CONTINUE
      DO 20 L=1,N 
         DR(L) =-TEMP*(SUM + FLOAT(L))
20    CONTINUE
                                  RETURN
      END 
C******************* WAT 2 ***********************************
C 
      SUBROUTINE F5(N,V,X) 
      REAL V(N),X(N+1)
      INTEGER N 
      COMMON KJAC,KFUNC 
C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION 
C     WE WISH TO CONTINUE ALONG. X CONTAINS THE VALUE OF (X(N),R) 
C     AT WHICH F IS COMPUTED.  V(N)  CONTAINS F(X,R) ON RETURN. 
C     X,N ARE UNCHANGED.
C 
C******************************************************************** 
C
C      SUBROUTINE F(X,V) -- EVALUATES WAT 2 FUNCTION AT THE POINT 
C         X, AND RETURNS THE VALUE IN V.
C
C******************************************************************** 
C
      KFUNC = KFUNC + 1 
      V(1) = 0.01*(X(1)+X(2)+1.)**3 
      DO 10 L=2,N-1 
         V(L) = 0.01*(X(L-1)+X(L)+X(L+1)+1.)**3 
10    CONTINUE
      V(N) = 0.01*(X(N-1)+X(N)+1.)**3 
      XLMDA = X(N+1)
      DO 20 L=1,N 
         V(L) = X(L) - XLMDA*V(L) 
20    CONTINUE
                       RETURN 
      END 
C* 
      SUBROUTINE DZF5(N,MXDM,DZ,X) 
      REAL DZ(MXDM,N), X(N+1) 
      INTEGER N,MXDM
      COMMON KJAC,KFUNC 
C ***  COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F
C      WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED
C      IN X(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N).
C 
C 
      KJAC = KJAC+1 
C 
C******************************************************************** 
C
C         THE JACOBIAN MATRIX FOR WAT 2 FUNCTION EVALUATED AT 
C
C******************************************************************** 
C
      DO 20 I=1,N 
      DO 10 L=1,N 
       DZ(L,I)= 0.0 
10    CONTINUE
20     CONTINUE 
       DZ(1,1)=0.03*(X(1)+X(2)+1.)**2
       DZ(1,2)=0.03*(X(1)+X(2)+1.)**2 
      DZ(N,N-1)=0.03*(X(N-1)+X(N)+1.)**2
      DZ(N,N)=0.03*(X(N-1)+X(N)+1.)**2
      DO 30 K=2,N-1 
        TEMP=0.03*(X(K-1)+X(K+1)+X(K)+1.)**2
         DZ(K,K-1) = TEMP 
        DZ(K,K)=TEMP
        DZ(K+1,K)=TEMP
30      CONTINUE
      XLMDA = X(N+1)
      DO 40 L=1,N 
         DO 35 J=1,N
            DZ(L,J) = -XLMDA*DZ(L,J)
35       CONTINUE 
      DZ(L,L) = 1.+DZ(L,L)
40    CONTINUE
       RETURN
       END
C 
      SUBROUTINE DRF5(N,DR,X)
      REAL DR(N),X(N+1) 
      INTEGER N 
C ***  COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F
C      WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN X(N+1)
C      THE DERIVATIVE IS RETURNED IN DR(N). 
      DR(1) = -0.01*(X(1)+X(2)+1.)**3 
      DO 10 L=2,N-1 
         DR(L) = -0.01*(X(L-1)+X(L)+X(L+1)+1.)**3 
10    CONTINUE
      DR(N) = -0.01*(X(N-1)+X(N)+1.)**3 
                       RETURN 
      END 
C ***************** BR ********************** 
      SUBROUTINE F6(N,V,X) 
      REAL V(N),X(N+1)
      INTEGER N 
      COMMON KJAC,KFUNC 
C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION 
C     WE WISH TO CONTINUE ALONG. X CONTAINS THE VALUE OF (X(N),R) 
C     AT WHICH F IS COMPUTED.  V(N)  CONTAINS F(X,R) ON RETURN. 
C     X,N ARE UNCHANGED.
C 
C******************************************************************** 
C
C      SUBROUTINE F(X,V) -- EVALUATES BROWN'S FUNCTION AT THE POINT
C         X, AND RETURNS THE VALUE IN V.
C
C******************************************************************** 
C
      KFUNC = KFUNC + 1 
       PROD=1.0E0 
       DO 10 J=1,N
10     PROD=PROD*X(J)
       V(1)=PROD-1.0E0
       SUM=0.0E0
       DO 20 J=1,N
20     SUM=SUM+X(J) 
       SUM=SUM-FLOAT(N+1) 
       DO 30 J=2,N
30     V(J)=SUM+X(J)
      TEMP = X(N+1) 
      DO 40 L=1,N 
         V(L)  = TEMP*V(L)+(1.-TEMP)*X(L) 
40    CONTINUE
       RETURN
       END
C 
C 
      SUBROUTINE DZF6(N,MXDM,DZ,X) 
      REAL DZ(MXDM,N), X(N+1) 
      INTEGER N,MXDM
      COMMON KJAC,KFUNC 
C ***  COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F
C      WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED
C      IN X(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N).
C 
C 
      KJAC = KJAC+1 
C 
C******************************************************************** 
C
C         THE JACOBIAN MATRIX FOR BROWN'S FUNCTION EVALUATED AT
C         THE POINT X, RETURNING THE VALUE IN V.
C
C******************************************************************** 
C
      DO 30 K=1,N 
       PROD=1.0E0 
       DO 10 J=1,K-1
10     PROD=PROD*X(J)
       DO 15 J=K+1,N
15     PROD=PROD*X(J)
      DZ(1,K)=PROD
       DO 20 J=2,N
20     DZ(J,K)=1.0
       IF (K .GT. 1) DZ(K,K)=DZ(K,K)+1.0
30    CONTINUE
      TEMP = X(N+1) 
      DO 50 I=1,N 
         DO 40 J=1,N
            DZ(J,I) = TEMP*DZ(J,I)
40       CONTINUE 
         DZ(I,I) = 1.-TEMP +DZ(I,I) 
50    CONTINUE
       RETURN
       END
      SUBROUTINE DRF6(N,DR,X)
      REAL DR(N),X(N+1) 
      INTEGER N 
C ***  COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F
C      WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN X(N+1)
C      THE DERIDRATIDRE IS RETURNED IN DR(N). 
       PROD=1.0E0 
       DO 10 J=1,N
10     PROD=PROD*X(J)
       DR(1)=PROD-1.0E0 
       SUM=0.0E0
       DO 20 J=1,N
20     SUM=SUM+X(J) 
       SUM=SUM-FLOAT(N+1) 
       DO 30 J=2,N
30     DR(J)=SUM+X(J) 
      DO 40 J=1,N 
         DR(J) = DR(J)- X(J)
40    CONTINUE
       RETURN
       END
C *********** FREUDENSTEIN AND ROTH -- REGULARIZING  *******
      SUBROUTINE F7(N,RESULT,ARG)
C ** FREUDENSTIEN AND ROTH ** 
      REAL RESULT(N),ARG(N+1) 
      INTEGER N 
      COMMON KJAC,KFUNC 
      KFUNC = KFUNC+1 
C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION 
C     WE WISH TO CONTINUE ALONG. ARG CONTAINS THE VALUE OF (X(N),R) 
C     AT WHICH F IS COMPUTED.  RESULT(N)  CONTAINS F(X,R) ON RETURN.
C     ARG,N ARE UNCHANGED.
C 
C 
      RESULT(1)=ARG(1)+5.*ARG(2)**2-ARG(2)**3-2.*ARG(2)-13. 
      RESULT(2)=ARG(1)+ARG(2)**2+ARG(2)**3-14.*ARG(2)-29. 
      XL = ARG(3) 
      RESULT(1) = XL*RESULT(1) + (1.-XL)*(ARG(1) - 15.) 
      RESULT(2) = XL*RESULT(2) + (1.-XL)*(ARG(2) +2.) 
                                        RETURN
      END 
C 
C 
      SUBROUTINE DZF7(N,MXDM,DZ,ARG) 
      REAL DZ(MXDM,N), ARG(N+1) 
      INTEGER N,MXDM
      COMMON KJAC,KFUNC 
      KJAC = KJAC + 1 
C ***  COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F
C      WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED
C      IN ARG(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N).
C 
      DZ(1,1) =  1.0
      DZ(2,1) = 1.0 
      DZ(1,2) = -3.*ARG(2)**2+10.*ARG(2)-2. 
      DZ(2,2) = 3.*ARG(2)**2+2.*ARG(2)-14.
      DO 20 J=1,2 
         DO 10 I=1,2
            DZ(I,J) = ARG(3)*DZ(I,J)
10       CONTINUE 
         DZ(J,J) = DZ(J,J) + 1.-ARG(3)
20    CONTINUE
                                       RETURN 
      END 
C 
C 
C 
      SUBROUTINE DRF7(N,DR,ARG)
      REAL DR(N),ARG(N+1) 
      INTEGER N 
C ***  COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F
C      WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN ARG(N+1)
C      THE DERIVATIVE IS RETURNED IN DR(N). 
C 
      DR(1)=ARG(1)+5.*ARG(2)**2-ARG(2)**3-2.*ARG(2)-13. 
      DR(2)=ARG(1)+ARG(2)**2+ARG(2)**3-14.*ARG(2)-29. 
      DR(1) = DR(1) - ARG(1) + 15.
      DR(2) = DR(2) - ARG(2) - 2. 
                                       RETURN 
      END
C ************ WATSON CURVE *******************
C
      SUBROUTINE FF(N,RESULT,ARG) 
      REAL RESULT(N),ARG(N+1) 
      INTEGER N 
C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION 
C     WE WISH TO CONTINUE ALONG. ARG CONTAINS THE VALUE OF (X(N),R) 
C     AT WHICH F IS COMPUTED.  RESULT(N)  CONTAINS F(X,R) ON RETURN.
C     ARG,N ARE UNCHANGED.
C 
C 
      SUM = 0.0 
      DO 10 L=1,N 
         SUM = SUM + ARG(L) 
10    CONTINUE
      DO 20 L=1,N 
      RESULT(L)=EXP(COS(L*SUM)) 
20    CONTINUE
                                        RETURN
      END 
C 
      SUBROUTINE F8(N,RESULT,ARG)
      REAL RESULT(N),ARG(N+1) 
      INTEGER N 
      COMMON KJAC,KFUNC 
      KFUNC = KFUNC + 1 
C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION 
C     WE WISH TO CONTINUE ALONG. ARG CONTAINS THE VALUE OF (X(N),R) 
C     AT WHICH F IS COMPUTED.  RESULT(N)  CONTAINS F(X,R) ON RETURN.
C     ARG,N ARE UNCHANGED.
C 
C 
       CALL FF(N,RESULT,ARG)
      DO 10 L=1,N 
         RESULT(L) = ARG(L) - ARG(N+1)*RESULT(L)
10    CONTINUE
                                      RETURN
      END 
C 
C 
C 
      SUBROUTINE DZF8(N,MXDM,DZ,ARG) 
      REAL DZ(MXDM,N), ARG(N+1),FN(1001)
      INTEGER N,MXDM
      COMMON KJAC,KFUNC 
      KJAC = KJAC+1 
C ***  COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F
C      WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED
C      IN ARG(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N).
C 
      SUM = 0.
      DO 10 L=1,N 
         SUM = SUM + ARG(L) 
10    CONTINUE
      CALL FF(N,FN,ARG) 
      DO 20 I=1,N 
         TEMP = I*FN(I)*SIN(I*SUM)*ARG(N+1) 
         DO 15 J=1,N
            DZ(I,J) =TEMP 
15       CONTINUE 
         DZ(I,I) = DZ(I,I) + 1. 
20    CONTINUE
                                       RETURN 
      END 
C 
C 
C 
      SUBROUTINE DRF8(N,DR,ARG)
      REAL DR(N),ARG(N+1) 
      INTEGER N 
C ***  COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F
C      WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN ARG(N+1)
C      THE DERIVATIVE IS RETURNED IN DR(N). 
C 
      CALL FF(N,DR,ARG) 
      DO 10 I=1,N 
         DR(I) = -DR(I) 
10    CONTINUE
                                       RETURN 
      END 
      SUBROUTINE F9(N,RESULT,ARG)
      PARAMETER(A=-25.,B=-25.)
      REAL RESULT(N),ARG(N+1) 
      INTEGER N 
      COMMON KJAC,KFUNC 
C *** THIS SUBROUTINE COMPUTES F(X,R) WHERE F IS THE FUNCTION 
C     WE WISH TO CONTINUE ALONG. ARG CONTAINS THE VALUE OF (X(N),R) 
C     AT WHICH F IS COMPUTED.  RESULT(N)  CONTAINS F(X,R) ON RETURN.
C     ARG,N ARE UNCHANGED.
C 
C 
C 
      KFUNC = KFUNC + 1 
      RESULT(1) = ARG(1)*EXP(A*ARG(3))
      RESULT(2) = ARG(2)*EXP(B*ARG(3))
                  RETURN
      END 
C 
C 
C 
      SUBROUTINE DZF9(N,MXDM,DZ,ARG) 
      PARAMETER(A=-25.,B=-25.)
      REAL DZ(MXDM,N), ARG(N+1) 
      INTEGER N,MXDM
      COMMON KJAC,KFUNC 
C ***  COMPUTES THE DERIVATIVE(NXN MATRIX) OF THE FUNCTION F
C      WITH RESPECT TO U AT THE VECTOR (U(N),R), WHICH IS PASSED
C      IN ARG(N+1). THE DERIVATIVE IS RETURNED IN DZ(N,N).
C 
C 
      KJAC = KJAC+1 
C 
      DO 20 J=1,N 
         DO 10 I=1,N
            DZ(I,J) = 0.
10       CONTINUE 
20    CONTINUE
C 
      DZ(1,1) =EXP(A*ARG(3))
      DZ(2,2) = EXP(B*ARG(3)) 
                                     RETURN 
      END 
C 
C 
C 
      SUBROUTINE DRF9(N,DR,ARG)
      PARAMETER(A=-25.,B=-25.)
      REAL DR(N),ARG(N+1) 
      INTEGER N 
C ***  COMPUTES THE DERIVATIVE(N-VECTOR) OF THE FUNCTION F
C      WITH RESPECT TO R AT (U,R) WHICH IS CONTAINED IN ARG(N+1)
C      THE DERIVATIVE IS RETURNED IN DR(N). 
C 
      DR(1) = ARG(1)*A*EXP(A*ARG(3))
      DR(2) = B*ARG(2)*EXP(B*ARG(3))
                                  RETURN
      END