C      ALGORITHM 688, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 17, NO. 2, PP. 153-166.  JUNE, 1991.


      SUBROUTINE PDECOL(T0,TOUT,DT,XBKPT,EPS,NINT,KORD,NCC,NPDE,MF,
     *                  INDEX,WORK,IWORK)
      IMPLICIT REAL*8(A-H,O-Z)
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C THIS IS THE AUG 1987 VERSION OF PDECOL, WHICH MAKES USE OF THE ALMOST
C BLOCK-DIAGONAL LINEAR EQUATION SOLVER COLROW (COMPRISING CRDCMP, FOR
C THE FACTORIZATION, AND CRSLVE FOR THE SOLUTION PHASE).
C
C THIS PACKAGE WAS CONSTRUCTED SO AS TO CONFORM TO AS MANY ANSI-FORTRAN
C RULES AS WAS CONVIENTLY POSSIBLE.  THE FORTRAN USED VIOLATES ANSI
C STANDARDS IN THE TWO WAYS LISTED BELOW....
C
C   1. SUBSCRIPTS OF THE GENERAL FORM C*V1 + V2 + V3 ARE USED
C      (POSSIBLY IN A PERMUTED ORDER), WHERE C IS AN INTEGER CONSTANT
C      AND V1, V2, AND V3 ARE INTEGER VARIABLES.
C
C   2. ARRAY NAMES APPEAR SINGLY IN DATA STATEMENTS IN THE ROUTINES
C      BSPLVN AND COSET.
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
C PDECOL IS THE DRIVER ROUTINE FOR A SOPHISTICATED PACKAGE OF
C SUBROUTINES WHICH IS DESIGNED TO SOLVE THE GENERAL SYSTEM OF
C NPDE NONLINEAR PARTIAL DIFFERENTIAL EQUATIONS OF AT MOST SECOND
C ORDER ON THE INTERVAL (XLEFT,XRIGHT) FOR T .GT. T0 WHICH IS OF THE
C FORM....
C
C      DU/DT  =  F( T, X, U, UX, UXX )
C
C WHERE
C
C          U  =  (  U(1),  U(2), ... ,  U(NPDE) )
C         UX  =  ( UX(1), UX(2), ... , UX(NPDE) )
C        UXX  =  (UXX(1),UXX(2), ... ,UXX(NPDE) ) .
C
C EACH U(K) IS A FUNCTION OF THE SCALAR QUANTITIES T AND X.
C UX(K) REPRESENTS THE FIRST PARTIAL DERIVATIVE OF U(K) WITH RESPECT
C TO THE VARIABLE X,  UXX(K) REPRESENTS THE SECOND PARTIAL DERIVATIVE
C OF U(K) WITH RESPECT TO THE VARIABLE X, AND DU/DT IS THE VECTOR OF
C PARTIAL DERIVATIVES OF U WITH RESPECT TO THE TIME VARIABLE T.
C F  REPRESENTS AN ARBITRARY VECTOR VALUED FUNCTION WHOSE NPDE
C COMPONENTS DEFINE THE RESPECTIVE PARTIAL DIFFERENTIAL EQUATIONS OF
C THE PDE SYSTEM.  SEE SUBROUTINE F DESCRIPTION BELOW.
C
C BOUNDARY CONDITIONS
C
C   IN THE NEW VERSION OF PDECOL, 2 BOUNDARY CONDITIONS
C   ARE REQUIRED FOR EACH PDE IN THE SYSTEM. THESE ARE IMPOSED AT XLEFT
C   AND RIGHT AND EACH MUST BE OF THE FORM....
C
C        B(U,UX)  =  Z(T)
C
C   WHERE  B  AND  Z  ARE ARBITRARY VECTOR VALUED FUNCTIONS WITH
C   NPDE COMPONENTS AND  U, UX, AND T  ARE AS ABOVE.  THESE BOUNDARY
C   CONDITIONS MUST BE CONSISTENT WITH THE INITIAL CONDITIONS WHICH ARE
C   DESCRIBED NEXT.
C
C INITIAL CONDITIONS
C
C   EACH SOLUTION COMPONENT  U(K)  IS ASSUMED TO BE A KNOWN (USER
C   PROVIDED) FUNCTION OF  X  AT THE INITIAL TIME T = T0.  THE
C   INITIAL CONDITION FUNCTIONS MUST BE CONSISTENT WITH THE BOUNDARY
C   CONDITIONS ABOVE, I.E. THE INITIAL CONDITION FUNCTIONS MUST
C   SATISFY THE BOUNDARY CONDITIONS FOR T = T0.  SEE SUBROUTINE UINIT
C   DESCRIPTION BELOW.
C-----------------------------------------------------------------------
C
C REQUIRED USER SUPPLIED SUBROUTINES
C
C THE USER IS REQUIRED TO CONSTRUCT THREE SUBPROGRAMS AND A MAIN
C PROGRAM WHICH DEFINE THE PDE PROBLEM WHOSE SOLUTION IS TO BE
C ATTEMPTED.  THE THREE SUBPROGRAMS ARE...
C
C 1)  SUBROUTINE F( T, X, U, UX, UXX, FVAL, NPDE )
C     DIMENSION U(NPDE), UX(NPDE), UXX(NPDE), FVAL(NPDE)
C        THIS ROUTINE DEFINES THE DESIRED PARTIAL DIFFERENTIAL
C        EQUATIONS TO BE SOLVED.  THE PACKAGE PROVIDES VALUES OF THE
C        INPUT SCALARS T AND X AND INPUT ARRAYS (LENGTH NPDE) U, UX,
C        AND UXX, AND THE USER MUST CONSTRUCT THIS ROUTINE TO COMPUTE
C        THE OUTPUT ARRAY FVAL (LENGTH NPDE) WHICH CONTAINS THE
C        CORRESPONDING VALUES OF THE RIGHT HAND SIDES OF THE DESIRED
C        PARTIAL DIFFERENTIAL EQUATIONS, I.E.
C
C        FVAL(K) = THE VALUE OF THE RIGHT HAND SIDE OF THE K-TH PDE IN
C                  THE PDE SYSTEM ABOVE, FOR K = 1 TO NPDE.
C
C        THE INCOMING VALUE OF THE SCALAR QUANTITY X WILL BE A
C        COLLOCATION POINT VALUE (SEE INITAL AND COLPNT) AND THE
C        INCOMING VALUES IN THE ARRAYS U, UX AND UXX CORRESPOND TO THIS
C        POINT X AND TIME T.
C     RETURN
C     END
C
C 2)  SUBROUTINE BNDRY( T, X, U, UX, DBDU, DBDUX, DZDT, NPDE )
C     DIMENSION U(NPDE), UX(NPDE), DZDT(NPDE)
C     DIMENSION DBDU(NPDE,NPDE), DBDUX(NPDE,NPDE)
C        THIS ROUTINE IS USED TO PROVIDE THE PDE PACKAGE WITH NEEDED
C        INFORMATION ABOUT THE BOUNDARY CONDITION FUNCTIONS B AND Z
C        ABOVE.  THE PACKAGE PROVIDES VALUES OF THE INPUT VARIABLES
C        T, X, U, AND UX, AND THE USER IS TO DEFINE THE CORRESPONDING
C        OUTPUT VALUES OF THE DERIVATIVES OF THE FUNCTIONS B AND Z
C        WHERE....
C           DBDU(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE
C                       VECTOR FUNCTION B(U,UX) ABOVE WITH RESPECT TO
C                       THE J-TH VARIABLE U(J).
C          DBDUX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE
C                       VECTOR FUNCTION B(U,UX) ABOVE WITH RESPECT TO
C                       THE J-TH VARIABLE UX(J).
C             DZDT(K) = DERIVATIVE OF THE K-TH COMPONENT OF THE VECTOR
C                       FUNCTION Z(T) ABOVE WITH RESPECT TO THE
C                       VARIABLE T.
C        NOTE... THE INCOMING VALUE OF X WILL BE EITHER XLEFT OR XRIGHT.
C        WE REQUIRE A BOUNDARY CONDITION TO BE IMPOSED AT BOTH ENDPOINTS
C        OF EACH PDE IN THE SYSTEM. THUS, FOR SAY THE K-TH PDE, WE REQUI
C        AT BOTH OF THE ENDPOINTS, XLEFT AND XRIGHT, THAT ONE OF DBDU(K,
C        OR DBDUX(K,K) SHOULD BE NON-ZERO WHEN BNDRY IS CALLED.
C        THIS ROUTINE CAN BE STRUCTURED AS FOLLOWS...
C        THE COMMON BLOCK /ENDPT/ IS NOT A PART OF PDECOL AND
C        MUST BE SUPPLIED AND DEFINED BY THE USER.
C     COMMON /ENDPT/ XLEFT
C     IF( X .NE. XLEFT ) GO TO 10
C        HERE DEFINE AND SET PROPER VALUES FOR DBDU(K,J), DBDUX(K,J),
C        AND DZDT(K) FOR K,J = 1 TO NPDE FOR THE LEFT BOUNDARY POINT
C        X = XLEFT.
C     RETURN
C  10 CONTINUE
C        HERE DEFINE AND SET PROPER VALUES FOR DBDU(K,J), DBDUX(K,J),
C        AND DZDT(K) FOR K,J = 1 TO NPDE FOR THE RIGHT BOUNDARY POINT
C        X = XRIGHT.
C     RETURN
C     END
C
C 3)  SUBROUTINE UINIT( X, U, NPDE )
C     DIMENSION U(NPDE)
C        THIS ROUTINE IS USED TO PROVIDE THE PDE PACKAGE WITH THE
C        NEEDED INITIAL CONDITION FUNCTION VALUES.  THE PACKAGE
C        PROVIDES A VALUE OF THE INPUT VARIABLE X, AND THE USER IS TO
C        DEFINE THE PROPER INITIAL VALUES (AT T = T0) FOR ALL OF THE
C        PDE COMPONENTS, I.E.
C           U(K) = DESIRED INITIAL VALUE OF PDE COMPONENT U(K) AT
C                  X AND T = T0 FOR K = 1 TO NPDE.
C        NOTE... THE INCOMING VALUE OF X WILL BE A COLLOCATION POINT
C        VALUE.  THE INITIAL CONDITIONS AND BOUNDARY CONDITIONS
C        MUST BE CONSISTENT (SEE ABOVE).
C     RETURN
C     END
C-----------------------------------------------------------------------
C
C OPTIONAL USER SUPPLIED SUBROUTINE
C
C IF THE USER DESIRES TO USE THE MF = 11 OR 21 OPTION IN ORDER TO SAVE
C ABOUT 10-20 PERCENT IN EXECUTION TIME (SEE BELOW), THEN THE USER MUST
C PROVIDE THE FOLLOWING SUBROUTINE WHICH PROVIDES INFORMATION ABOUT THE
C DERIVATIVES OF THE FUNCTION F ABOVE. THIS PROVIDES FOR MORE EFFICIENT
C JACOBIAN MATRIX GENERATION.  ON MOST COMPUTER SYSTEMS, THE USER WILL
C BE REQUIRED TO SUPPLY THIS SUBROUTINE AS A DUMMY SUBROUTINE IF THE
C OPTIONS MF = 12 OR 22 ARE USED (SEE BELOW).
C
C 1)  SUBROUTINE DERIVF( T, X, U, UX, UXX, DFDU, DFDUX, DFDUXX, NPDE )
C     DIMENSION U(NPDE), UX(NPDE), UXX(NPDE)
C     DIMENSION DFDU(NPDE,NPDE), DFDUX(NPDE,NPDE), DFDUXX(NPDE,NPDE)
C        THE PACKAGE PROVIDES VALUES OF THE INPUT VARIABLES T, X, U, UX,
C        AND UXX, AND THE USER SHOULD CONSTRUCT THIS ROUTINE TO PROVIDE
C        THE FOLLOWING CORRESPONDING VALUES OF THE OUTPUT ARRAYS
C        DFDU, DFDUX, AND DFDUXX FOR K,J = 1 TO NPDE...
C           DFDU(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE
C                       PDE DEFINING FUNCTION  F  WITH RESPECT TO THE
C                       VARIABLE U(J).
C          DFDUX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE
C                       PDE DEFINING FUNCTION  F  WITH RESPECT TO THE
C                       VARIABLE UX(J).
C         DFDUXX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE
C                       PDE DEFINING FUNCTION  F  WITH RESPECT TO THE
C                       VARIABLE UXX(J).
C        NOTE... THE INCOMING VALUE OF  X  WILL BE A COLLOCATION POINT
C        VALUE.
C     RETURN
C     END
C-----------------------------------------------------------------------
C
C METHODS USED
C
C   THE PACKAGE PDECOL IS BASED ON THE METHOD OF LINES AND USES A
C   FINITE ELEMENT COLLOCATION PROCEDURE (WITH PIECEWISE POLYNOMIALS
C   AS THE TRIAL SPACE) FOR THE DISCRETIZATION OF THE SPATIAL VARIABLE
C   X.  THE COLLOCATION PROCEDURE REDUCES THE PDE SYSTEM TO A SEMI-
C   DISCRETE SYSTEM WHICH THEN DEPENDS ONLY ON THE TIME VARIABLE T.
C   THE TIME INTEGRATION IS THEN ACCOMPLISHED BY USE OF SLIGHTLY
C   MODIFIED STANDARD TECHNIQUES (SEE REFS. 1,2).
C
C   PIECEWISE POLYNOMIALS
C
C   THE USER IS REQUIRED TO SELECT THE PIECEWISE POLYNOMIAL SPACE
C   WHICH IS TO BE USED TO COMPUTE HIS APPROXIMATE SOLUTION. FIRST, THE
C   ORDER, KORD, OF THE POLYNOMIALS TO BE USED MUST BE SPECIFIED
C   (KORD = POLYNOMIAL DEGREE + 1).  NEXT, THE NUMBER OF PIECES
C   (INTERVALS), NINT, INTO WHICH THE SPATIAL DOMAIN (XLEFT,XRIGHT) IS
C   TO BE DIVIDED, IS CHOSEN.  THE NINT + 1 DISTINCT BREAKPOINTS OF
C   THE DOMAIN MUST BE DEFINED AND SET INTO THE ARRAY XBKPT IN
C   STRICTLY INCREASING ORDER, I.E.
C   XLEFT=XBKPT(1) .LT. XBKPT(2) .LT. ... .LT. XBKPT(NINT+1)=XRIGHT.
C   THE APPROXIMATE SOLUTION AT ANY TIME T WILL BE A POLYNOMIAL OF
C   ORDER KORD OVER EACH SUBINTERVAL (XBKPT(I),XBKPT(I+1)).  THE
C   NUMBER OF CONTINUITY CONDITIONS, NCC, TO BE IMPOSED ACROSS ALL OF
C   THE BREAKPOINTS IS THE LAST PIECE OF USER SUPPLIED DATA WHICH IS
C   REQUIRED TO UNIQUELY DETERMINE THE DESIRED PIECEWISE POLYNOMIAL
C   SPACE.  FOR EXAMPLE, NCC = 2 WOULD REQUIRE THAT THE APPROXIMATE
C   SOLUTION (MADE UP OF THE SEPARATE POLYNOMIAL PIECES) AND ITS FIRST
C   SPATIAL DERIVATIVE BE CONTINUOUS AT THE BREAKPOINTS AND HENCE ON
C   THE ENTIRE DOMAIN (XLEFT,XRIGHT).  NCC = 3 WOULD REQUIRE THAT THE
C   APPROXIMATE SOLUTION AND ITS FIRST AND SECOND SPATIAL DERIVATIVES
C   BE CONTINUOUS AT THE BREAKPOINTS, ETC. THE DIMENSION OF THIS LINEAR
C   SPACE IS KNOWN AND FINITE AND IS NCPTS = KORD*NINT - NCC*(NINT-1).
C   THE WELL-KNOWN B-SPLINE BASIS (SEE REF. 3) FOR THIS SPACE IS USED
C   BY PDECOL AND IT CONSISTS OF NCPTS KNOWN PIECEWISE POLYNOMIAL
C   FUNCTIONS BF(I,X), FOR I=1 TO NCPTS, WHICH DO NOT DEPEND ON THE
C   TIME VARIABLE T. WE WISH TO EMPHASIZE THAT THE PIECEWISE POLYNOMIAL
C   SPACE USED IN PDECOL (WHICH IS SELECTED BY THE USER) WILL DETERMINE
C   THE MAGNITUDE OF THE SPATIAL DISCRETIZATION ERRORS IN THE COMPUTED
C   APPROXIMATE SOLUTION.  THE PACKAGE HAS NO CONTROL OVER ERRORS
C   INTRODUCED BY THE USERS CHOICE OF THIS SPACE.  SEE INPUT PARAMETERS
C   BELOW.
C
C   COLLOCATION OVER PIECEWISE POLYNOMIALS
C
C   THE BASIC ASSUMPTION MADE IS THAT THE APPROXIMATE SOLUTION
C   SATISFIES
C                       NCPTS
C             U(T,X)  =  SUM  C(I,T) * BF(I,X)
C                        I=1
C
C   WHERE THE UNKNOWN COEFFICIENTS C DEPEND ONLY ON THE TIME T AND
C   THE KNOWN BASIS FUNCTIONS DEPEND ONLY ON X (WE HAVE ASSUMED THAT
C   NPDE = 1 FOR CONVENIENCE).  SO, AT ANY GIVEN TIME T THE APPROX-
C   IMATE SOLUTION IS A PIECEWISE POLYNOMIAL IN THE USER CHOSEN SPACE.
C   THE SEMI-DISCRETE EQUATIONS (ACTUALLY ORDINARY DIFFERENTIAL
C   EQUATIONS) WHICH DETERMINE THE COEFFICIENTS C ARE OBTAINED BY
C   REQUIRING THAT THE ABOVE APPROXIMATE U(T,X) SATISFY THE PDE AND
C   BOUNDARY CONDITIONS EXACTLY AT A SET OF NCPTS COLLOCATION POINTS
C   (SEE COLPNT).  THUS, PDECOL ACTUALLY COMPUTES THE BASIS FUNCTION
C   COEFFICIENTS RATHER THAN SPECIFIC APPROXIMATE SOLUTION VALUES.
C
C-----------------------------------------------------------------------
C
C   THE MODIFICATION MADE TO THE ORIGINAL PDECOL PROGRAM CONSISTS
C   ENTIRELY OF THE REMOVAL OF THE BAND LINEAR EQUATION SOLVER, DECB
C   (FOR THE DECOMPOSITION PHASE)  AND SOLB (FOR THE SOLUTION PHASE),
C   AND ITS SUBSTITUTION BY CRDCMP, CRSLVE, THE CONSTITUENT PARTS OF
C   COLROW [5]. THE CHANGES NECESSITATED IN SUBROUTINE PDECOL BY THIS
C   SUBSTITUTION ARE LIMITED TO :
C
C          (1)  CHANGES TO THE DIMENSION OF PW OR WORK(IW17).
C          (2)  RESULTING CHANGES IN THE CHECKS ON INPUT PARAMETERS.
C          (3)  THE REMOVAL OF SOME VARIABLES USED SOLELY TO PROVIDE
C               INFORMATION TO DECB, SOLB, AND THE INSERTION OF NEW
C               VARIABLES TO PASS INFORMATION TO CRDCMP, CRSLVE.
C
C-----------------------------------------------------------------------
C
C   REFERENCES
C
C   1. MADSEN, N.K. AND R.F. SINCOVEC, PDECOL - COLLOCATION SOFTWARE
C        FOR PARTIAL DIFFERENTIAL EQUATIONS, ACM-TOMS, VOL. 5, NO. 3,
C        1979, PP. 326-351.
C   2. SINCOVEC, R.F. AND N.K. MADSEN, SOFTWARE FOR NONLINEAR PARTIAL
C        DIFFERENTIAL EQUATIONS, ACM-TOMS, VOL. 1, NO. 3,
C        SEPTEMBER 1975, PP. 232-260.
C   3. HINDMARSH, A.C., PRELIMINARY DOCUMENTATION OF GEARIB.. SOLUTION
C        OF IMPLICIT SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS WITH
C        BANDED JACOBIANS, LAWRENCE LIVERMORE LAB, UCID-30130, FEBRUARY
C        1976.
C   4. DEBOOR, C., PACKAGE FOR CALCULATING WITH B-SPLINES, SIAM J.
C        NUMER. ANAL., VOL. 14, NO. 3, JUNE 1977, PP. 441-472.
C   5. DIAZ, J.C., G. FAIRWEATHER, AND P.  KEAST, COLROW AND ARCECO :
C        FORTRAN PACKAGES FOR SOLVING CERTAIN ALMOST BLOCK DIAGONAL
C        LINEAR SYSTEMS BY MODIFIED ALTERNATE ROW AND COLUMN
C        ELIMINATION.  ACM TRAN. MATH. SOFTW. VOL. 9, NO. 3, SEPT. 1983,
C        PP. 376,380.
C
C-----------------------------------------------------------------------
C
C USE OF PDECOL
C
C PDECOL IS CALLED ONCE FOR EACH DESIRED OUTPUT VALUE (TOUT) OF THE
C TIME T, AND IT IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR,
C STIFIB, WHICH ADVANCES THE TIME BY TAKING SINGLE STEPS UNTIL
C T .GE. TOUT.  INTERPOLATION TO THE EXACT TIME TOUT IS THEN DONE.
C SEE TOUT BELOW.
C
C
C SUMMARY OF SUGGESTED INPUT VALUES
C
C   IT IS OF COURSE IMPOSSIBLE TO SUGGEST INPUT PARAMETER VALUES WHICH
C   ARE APPROPRIATE FOR ALL PROBLEMS.  THE FOLLOWING SUGGESTIONS ARE TO
C   BE USED ONLY IF YOU HAVE NO IDEA OF BETTER VALUES FOR YOUR PROBLEM.
C
C   DT    = 1.E-10
C   XBKPT = CHOOSE  NINT+1  EQUALLY SPACED VALUES SUCH THAT XBKPT(1) =
C           XLEFT AND  XBKPT(NINT+1) = XRIGHT.
C   EPS   = 1.E-4
C   NINT  = ENOUGH SO THAT ANY FINE STRUCTURE OF THE PROBLEM MAY BE
C           RESOLVED.
C   KORD  = 4
C   NCC   = 2
C   MF    = 22
C   INDEX = 1 (ON FIRST CALL ONLY, THEN 0 THEREAFTER).
C
C
C THE INPUT PARAMETERS ARE..
C   T0    =  THE INITIAL VALUE OF T, THE INDEPENDENT VARIABLE
C              (USED ONLY ON FIRST CALL).
C   TOUT  =  THE VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT.  SINCE
C              THE PACKAGE CHOOSES ITS OWN TIME STEP SIZES, THE
C              INTEGRATION WILL NORMALLY GO SLIGHTLY BEYOND TOUT
C              AND THE PACKAGE WILL INTERPOLATE TO T = TOUT.
C   DT    =  THE INITIAL STEP SIZE IN T, IF INDEX = 1, OR, THE
C              MAXIMUM STEP SIZE ALLOWED (MUST BE .GT. 0), IF INDEX = 3.
C              USED FOR INPUT ONLY WHEN INDEX = 1 OR 3. SEE BELOW.
C   XBKPT =  THE ARRAY OF PIECEWISE POLYNOMIAL BREAKPOINTS.
C              THE NINT+1 VALUES MUST BE STRICTLY INCREASING WITH
C              XBKPT(1) = XLEFT AND XBKPT(NINT+1) = XRIGHT (USED ONLY
C              ON FIRST CALL).
C   EPS   =  THE RELATIVE TIME ERROR BOUND  (USED ONLY ON THE
C              FIRST CALL, UNLESS INDEX = 4).  SINGLE STEP ERROR
C              ESTIMATES DIVIDED BY CMAX(I) WILL BE KEPT LESS THAN
C              EPS IN ROOT-MEAN-SQUARE NORM. THE VECTOR CMAX OF WEIGHTS
C              IS COMPUTED IN PDECOL.  INITIALLY CMAX(I) IS SET TO
C              ABS(C(I)), WITH A DEFAULT VALUE OF 1 IF ABS(C(I)) .LT. 1.
C              THEREAFTER, CMAX(I) IS THE LARGEST VALUE
C              OF ABS(C(I)) SEEN SO FAR, OR THE INITIAL CMAX(I) IF
C              THAT IS LARGER.  TO ALTER EITHER OF THESE, CHANGE THE
C              APPROPRIATE STATEMENTS IN THE DO-LOOPS ENDING AT
C              STATEMENTS 50 AND 130 BELOW.  THE USER SHOULD EXERCISE
C              SOME DISCRETION IN CHOOSING EPS.  IN GENERAL, THE
C              OVERALL RUNNING TIME FOR A PROBLEM WILL BE GREATER IF
C              EPS IS CHOSEN SMALLER. THERE IS USUALLY LITTLE REASON TO
C              CHOOSE EPS MUCH SMALLER THAN THE ERRORS WHICH ARE BEING
C              INTRODUCED BY THE USERS CHOICE OF THE POLYNOMIAL SPACE.
C              SEE RELATED COMMENTS CONCERNING CMAX BELOW STATEMENT 40.
C   NINT  =  THE NUMBER OF SUBINTERVALS INTO WHICH THE SPATIAL DOMAIN
C              (XLEFT,XRIGHT) IS TO BE DIVIDED (MUST BE .GE. 1)
C              (USED ONLY ON FIRST CALL).
C   KORD  =  THE ORDER OF THE PIECEWISE POLYNOMIAL SPACE TO BE USED.
C              ITS VALUE MUST BE GREATER THAN 2 AND LESS THAN 21.  FOR
C              FIRST ATTEMPTS WE RECOMMEND KORD = 4.  IF THE SOLUTION
C              IS SMOOTH AND MUCH ACCURACY IS DESIRED, HIGHER VALUES
C              MAY PROVE TO BE MORE EFFICIENT.  WE HAVE SELDOM USED
C              VALUES OF KORD IN EXCESS OF 8 OR 9, THOUGH THEY ARE
C              AVAILABLE FOR USE IN PDECOL (USED ONLY ON FIRST CALL).
C   NCC   =  THE NUMBER OF CONTINUITY CONDITIONS TO BE IMPOSED ON THE
C              APPROXIMATE SOLUTION AT THE BREAKPOINTS IN XBKPT.
C              NCC MUST BE GREATER THAN 1 AND LESS THAN KORD.  WE
C              RECOMMEND THE USE OF NCC = 2 (WITH NOGAUS = 0, SEE
C              BELOW), SINCE THEORY PREDICTS THAT DRAMATICALLY MORE
C              ACCURATE RESULTS CAN OFTEN BE OBTAINED USING THIS CHOICE
C              (USED ONLY ON FIRST CALL).
C   NPDE  =  THE NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS IN THE SYSTEM
C              TO BE SOLVED (USED ONLY ON FIRST CALL).
C   MF    =  THE METHOD FLAG  (USED ONLY ON FIRST CALL, UNLESS
C              INDEX = 4).  ALLOWED VALUES ARE 11, 12, 21, 22.
C              FOR FIRST ATTEMPTS WE RECOMMEND THE USE OF MF = 22.
C              MF HAS TWO DECIMAL DIGITS, METH AND MITER
C              (MF = 10*METH + MITER).
C              METH IS THE BASIC METHOD INDICATOR..
C                METH = 1  MEANS THE ADAMS METHODS (GENERALIZATIONS OF
C                          CRANK-NICOLSON).
C                METH = 2  MEANS THE BACKWARD DIFFERENTIATION
C                          FORMULAS (BDF), OR STIFF METHODS OF GEAR.
C              MITER IS THE ITERATION METHOD INDICATOR
C              AND DETERMINES HOW THE JACOBIAN MATRIX IS
C              TO BE COMPUTED..
C                MITER = 1 MEANS CHORD METHOD WITH ANALYTIC JACOBIAN.
C                          FOR THIS USER SUPPLIES SUBROUTINE DERIVF.
C                          SEE DESCRIPTION ABOVE.
C                MITER = 2 MEANS CHORD METHOD WITH JACOBIAN CALCULATED
C                          INTERNALLY BY FINITE DIFFERENCES.  SEE
C                          SUBROUTINES PSETIB AND DIFFF.
C   INDEX =  INTEGER USED ON INPUT TO INDICATE TYPE OF CALL,
C              WITH THE FOLLOWING VALUES AND MEANINGS..
C                 1    THIS IS THE FIRST CALL FOR THIS PROBLEM.
C                 0    THIS IS NOT THE FIRST CALL FOR THIS PROBLEM,
C                      AND INTEGRATION IS TO CONTINUE.
C                 2    SAME AS 0 EXCEPT THAT TOUT IS TO BE HIT
C                      EXACTLY (NO INTERPOLATION IS DONE).  SEE NOTE
C                      BELOW.  ASSUMES TOUT .GE. THE CURRENT T.
C                      IF TOUT IS .LT. THE CURRENT TIME, THEN TOUT IS
C                      RESET TO THE CURRENT TIME AND CONTROL IS
C                      RETURNED TO THE USER.  A CALL TO VALUES WILL
C                      PRODUCE ANSWERS FOR THE NEW VALUE OF TOUT.
C                 3    SAME AS 0 EXCEPT CONTROL RETURNS TO CALLING
C                      PROGRAM AFTER ONE STEP.  TOUT IS IGNORED AND
C                      DT MUST BE SET .GT. 0 TO A MAXIMUM ALLOWED
C                      DT VALUE. SEE ABOVE.
C                 4    THIS IS NOT THE FIRST CALL FOR THE PROBLEM,
C                      AND THE USER HAS RESET EPS AND/OR MF.
C              SINCE THE NORMAL OUTPUT VALUE OF INDEX IS 0,
C              IT NEED NOT BE RESET FOR NORMAL CONTINUATION.
C
C NOTE.. THE PACKAGE MUST HAVE TAKEN AT LEAST ONE SUCCESSFUL TIME
C STEP BEFORE A CALL WITH INDEX = 2 OR 4 IS ALLOWED.
C AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURRED AND A NORMAL
C CONTINUATION IS DESIRED, SIMPLY RESET TOUT AND CALL AGAIN.
C ALL OTHER PARAMETERS WILL BE READY FOR THE NEXT CALL.
C A CHANGE OF PARAMETERS WITH INDEX = 4 CAN BE MADE AFTER
C EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN PROVIDED AT LEAST
C ONE SUCCESSFUL TIME STEP HAS BEEN MADE.
C
C   WORK  =  FLOATING POINT WORKING ARRAY FOR PDECOL.  WE RECOMMEND
C          THAT IT BE INITIALIZED TO ZERO BEFORE THE FIRST CALL
C          TO PDECOL.  ITS TOTAL LENGTH MUST BE AT LEAST
C
C          KORD + 4*NPDE + NCPTS*(3*KORD + 2) +
C          NPDE*NCPTS*(MAXDER + 6) + (NPDE**2)*(13+KORD*(KORD-NCC)*NINT)
C
C              WHERE MAXDER IS DEFINED BELOW (SEE STORAGE
C              ALLOCATION).
C
C   IWORK =  INTEGER WORKING ARRAY FOR PDECOL.  THE FIRST TWO
C              LOCATIONS MUST BE DEFINED AS FOLLOWS...
C              IWORK(1) = LENGTH OF USERS ARRAY WORK
C              IWORK(2) = LENGTH OF USERS ARRAY IWORK
C              THE TOTAL LENGTH OF IWORK MUST BE AT LEAST
C              NCPTS*(NPDE + 1).
C OUTPUT
C
C THE SOLUTION VALUES ARE NOT RETURNED DIRECTLY TO THE USER BY PDECOL.
C THE METHODS USED IN PDECOL COMPUTE BASIS FUNCTION COEFFICIENTS, SO
C THE USER (AFTER A RETURN FROM PDECOL) MUST CALL THE PACKAGE ROUTINE
C VALUES TO OBTAIN HIS APPROXIMATE SOLUTION VALUES AT ANY DESIRED SPACE
C POINTS X AT THE TIME T = TOUT.  SEE THE COMMENTS IN SUBROUTINE VALUES
C FOR DETAILS ON HOW TO PROPERLY MAKE THE CALL.
C
C THE COMMON BLOCK /GEAR0/ CAN BE ACCESSED EXTERNALLY BY THE USER
C IF DESIRED.  IT CONTAINS THE STEP SIZE LAST USED (SUCCESSFULLY),
C THE ORDER LAST USED (SUCCESSFULLY), THE NUMBER OF STEPS TAKEN
C SO FAR, THE NUMBER OF RESIDUAL EVALUATIONS (RES CALLS) SO FAR,
C AND THE NUMBER OF MATRIX EVALUATIONS (PSETIB CALLS) SO FAR.
C DIFFUN CALLS ARE COUNTED IN WITH RESIDUAL EVALUATIONS.
C
C THE OUTPUT PARAMETERS ARE..
C   DT    =  THE STEP SIZE USED LAST, WHETHER SUCCESSFULLY OR NOT.
C   TOUT  =  THE OUTPUT VALUE OF T.  IF INTEGRATION WAS SUCCESSFUL,
C              AND THE INPUT VALUE OF INDEX WAS NOT 3, TOUT IS
C              UNCHANGED FROM ITS INPUT VALUE.  OTHERWISE, TOUT
C              IS THE CURRENT VALUE OF T TO WHICH THE INTEGRATION
C              HAS BEEN COMPLETED.
C   INDEX =  INTEGER USED ON OUTPUT TO INDICATE RESULTS,
C              WITH THE FOLLOWING VALUES AND MEANINGS..
C         0    INTEGRATION WAS COMPLETED TO TOUT OR BEYOND.
C        -1    THE INTEGRATION WAS HALTED AFTER FAILING TO PASS THE
C              ERROR TEST EVEN AFTER REDUCING DT BY A FACTOR OF
C              1.E10 FROM ITS INITIAL VALUE.
C        -2    AFTER SOME INITIAL SUCCESS, THE INTEGRATION WAS
C              HALTED EITHER BY REPEATED ERROR TEST FAILURES OR BY
C              A TEST ON EPS.  TOO MUCH ACCURACY HAS BEEN REQUESTED.
C        -3    THE INTEGRATION WAS HALTED AFTER FAILING TO ACHIEVE
C              CORRECTOR CONVERGENCE EVEN AFTER REDUCING DT BY A
C              FACTOR OF 1.E10 FROM ITS INITIAL VALUE.
C        -4    SINGULAR MATRIX ENCOUNTERED.  PROBABLY DUE TO STORAGE
C              OVERWRITES.
C        -5    INDEX WAS 4 ON INPUT, BUT THE DESIRED CHANGES OF
C              PARAMETERS WERE NOT IMPLEMENTED BECAUSE TOUT
C              WAS NOT BEYOND T.  INTERPOLATION TO T = TOUT WAS
C              PERFORMED AS ON A NORMAL RETURN.  TO TRY AGAIN,
C              SIMPLY CALL AGAIN WITH INDEX = 4 AND A NEW TOUT.
C        -6    ILLEGAL INDEX VALUE.
C        -7    ILLEGAL EPS VALUE.
C        -8    AN ATTEMPT TO INTEGRATE IN THE WRONG DIRECTION.  THE
C              SIGN OF DT IS WRONG RELATIVE TO T0 AND TOUT.
C        -9    DT .EQ. 0.0.
C       -10    ILLEGAL NINT VALUE.
C       -11    ILLEGAL KORD VALUE.
C       -12    ILLEGAL NCC VALUE.
C       -13    ILLEGAL NPDE VALUE.
C       -14    ILLEGAL MF VALUE.
C       -15    ILLEGAL BREAKPOINTS - NOT STRICTLY INCREASING.
C       -16    INSUFFICIENT STORAGE FOR WORK OR IWORK.
C-----------------------------------------------------------------------
C
C
C SUMMARY OF ALL PACKAGE ROUTINES
C
C PDECOL - STORAGE ALLOCATION, ERROR CHECKING, INITIALIZATION, REPEATED
C          CALLS TO STIFIB TO ADVANCE THE TIME.
C
C INTERP - INTERPOLATES COMPUTED BASIS FUNCTION COEFFICIENTS TO THE
C          DESIRED OUTPUT TIMES, TOUT, FOR USE BY VALUES.
C
C INITAL - INITIALIZATION, GENERATION AND STORAGE OF PIECEWISE POLY-
C          NOMIAL SPACE BASIS FUNCTION VALUES AND DERIVATIVES, DET-
C          ERMINES THE BASIS FUNCTION COEFFICINTS OF THE PIECEWISE
C          POLYNOMIALS WHICH INTERPOLATE THE USERS INITIAL CONDITIONS.
C
C COLPNT - GENERATION OF REQUIRED COLLOCATION POINTS.
C
C BSPLVD - B-SPLINE PACKAGE ROUTINES WHICH ALLOW FOR EVALUATION OF
C BSPLVN   ANY B-SPLINE BASIS FUNCTION OR DERIVATIVE VALUE.
C INTERV
C
C VALUES - GENERATION AT ANY POINT(S) OF VALUES OF THE COMPUTED
C          APPROXIMATE SOLUTION AND ITS DERIVATIVES WHICH ARE
C          PIECEWISE POLYNOMIALS.  THE SUBROUTINE IS CALLED ONLY BY
C          THE USER.
C
C STIFIB - CORE INTEGRATOR, TAKES SINGLE TIME STEPS TO ADVANCE THE
C          TIME.  ASSEMBLES AND SOLVES THE PROPER NONLINEAR EQUATIONS
C          WHICH ARE RELATED TO USE OF ADAMS OR GEAR TYPE INTEGRATION
C          FORMULAS.  CHOOSES PROPER STEP SIZE AND INTEGRATION FORMULA
C          ORDER TO MAINTAIN A DESIRED ACCURACY.  DESIGNED FOR ODE
C          PROBLEMS OF THE FORM A * (DY/DT) = G(T,Y).
C
C COSET  - GENERATES INTEGRATION FORMULA AND ERROR CONTROL COEFFICIENTS.
C
C RES    - COMPUTES RESIDUAL VECTORS USED IN SOLVING THE NONLINEAR
C          EQUATIONS BY A MODIFIED NEWTON METHOD.
C
C DIFFUN - COMPUTES A**-1 * G(T,Y) WHERE A AND G ARE AS ABOVE (STIFIB).
C
C ADDA   - ADDS THE A MATRIX TO A GIVEN MATRIX IN ALMOST BLOCK DIAGONAL
C
C EVAL   - EVALUATES THE COMPUTED PIECEWISE POLYNOMIAL SOLUTION AND
C          DERIVATIVES AT COLLOCATION POINTS.
C
C GFUN   - EVALUATES THE FUNCTION G(T,Y) BY CALLING EVAL AND THE USER
C          SUBROUTINES F AND BNDRY.
C
C PSETIB - GENERATES PROPER JACOBIAN MATRICES REQUIRED BY THE MODIFIED
C          NEWTON METHOD.
C
C DIFFF  - PERFORMS SAME ROLE AS THE USER ROUTINE DERIVF. COMPUTES
C          DERIVATIVE APPROXIMATIONS BY USE OF FINITE DIFFERENCES.
C
C CRDCMP - PERFORM AN LU DECOMPOSTION  OF THE ALMOST BLOCK DIAGONAL
C          LINEAR SYSTEMS ARISING.
C CRSLVE - PERFORM THE RIGHT-HAND SIDE MODIFICATIONS AND THE BACK-
C          SUBSTITUTIONS REQUIRED TO SOLVE THE ALMOST BLOCK DIAGONAL
C          SYSTEM AFTER IT HAS BEEN FACTORED.
C
C-----------------------------------------------------------------------
C
C
C STORAGE ALLOCATION
C
C SINCE PDECOL IS A DYNAMICALLY DIMENSIONED PROGRAM, MOST OF ITS
C WORKING STORAGE IS PROVIDED BY THE USER IN THE ARRAYS WORK AND IWORK.
C THE FOLLOWING GIVES A LIST OF THE ARRAYS WHICH MAKE UP THE CONTENTS
C WORK AND IWORK, THEIR LENGTHS, AND THEIR USES.  WHEN MORE THAN ONE
C NAME IS GIVEN, IT INDICATES THAT DIFFERENT NAMES ARE USED FOR THE
C SAME ARRAY IN DIFFERENT PARTS OF THE PROGRAM.  THE DIFFERENT NAMES
C OCCUR BECAUSE PDECOL IS AN AMALGAMATION OF SEVERAL OTHER CODES
C WRITTEN BY DIFFERENT PEOPLE AND WE HAVE TRIED TO LEAVE THE SEPARATE
C PARTS AS UNCHANGED FROM THEIR ORIGINAL VERSIONS AS POSSIBLE.
C
C
C   NAMES       LENGTH        USE
C ---------     ------------  -------------------------------------
C
C BC            4*NPDE**2     BOUNDARY CONDITION INFORMATION.
C WORK
C
C A             3*KORD*NCPTS  BASIS FUNCTION VALUES AT COLLOCATION POINT
C WORK(IW1)
C
C XT            NCPTS + KORD  BREAKPOINT SEQUENCE FOR GENERATION OF BASI
C WORK(IW2)                   FUNCTION VALUES.
C
C XC            NCPTS         COLLOCATION POINTS.
C WORK(IW3)
C
C CMAX          NPDE*NCPTS    VALUES USED IN ESTIMATING TIME
C YMAX                        INTEGRATION ERRORS.
C WORK(IW4)
C
C ERROR         NPDE*NCPTS    TIME INTEGRATION ERRORS.
C WORK(IW5)
C
C SAVE1         NPDE*NCPTS    WORKING STORAGE FOR THE TIME INTEGRATION
C WORK(IW6)                   METHOD.
C
C SAVE2         NPDE*NCPTS    WORKING STORAGE FOR THE TIME INTEGRATION
C WORK(IW7)                   METHOD.
C
C SAVE3         NPDE*NCPTS    WORKING STORAGE FOR THE TIME INTEGRATION
C WORK(IW8)                   METHOD.
C
C UVAL          3*NPDE        WORKING STORAGE FOR VALUES OF U, UX, AND
C WORK(IW9)                   USS AT ONE POINT.
C
C C             NPDE*NCPTS*   CURRENT BASIS FUNCTION COEFFICIENT VALUES
C Y               (MAXDER+1)  AND THEIR SCALED TIME DERIVATIVES.
C WORK(IW10)
C
C DFDU          NPDE**2       WORKING STORAGE USED TO COMPUTE THE
C WORK(IW11)                  JACOBIAN MATRIX.
C
C DFDUX         NPDE**2       WORKING STORAGE USED TO COMPUTE THE
C WORK(IW12)                  JACOBIAN MATRIX.
C
C DFDUXX        NPDE**2       WORKING STORAGE USED TO COMPUTE THE
C WORK(IW13)                  JACOBIAN MATRIX.
C
C DBDU          NPDE**2       BOUNDARY CONDITION INFORMATION.
C WORK(IW14)
C
C DBDUX         NPDE**2       BOUNDARY CONDITION INFORMATION.
C WORK(IW15)
C
C DZDT          NPDE          BOUNDARY CONDITION INFORMATION.
C WORK(IW16)
C
C PW            NPDE**2*(NCPTS*KORD-     STORAGE AND PROCESSING
C WORK(IW17)    NCC*(KORD-2))            OF THE JACOBIAN MATRIX.
C
C ILEFT         NCPTS         POINTERS TO BREAKPOINT SEQUENCE FOR
C IWORK                       GENERATION OF BASIS FUNCTION VALUES.
C
C IPIV          NPDE*NCPTS    PIVOT INFORMATION FOR THE LU DECOMPOSED
C IWORK(IW18)                 JACOBIAN MATRIX PW.
C
C WHERE...
C
C      NCPTS = KORD*NINT - NCC*(NINT-1)
C      MAXDER = 5 UNLESS OTHERWISE SET BY THE USER INTO /OPTION/.
C
C      IN THE ORIGINAL VERSION OF PDECOL A SWITCH CALLED IQUAD
C      (IQUAD = 1 IF KORD = 3 AND A NULL BOUNDARY CONDITION EXISTS
C      IQUAD = 0 OTHERWISE) WAS USED TO HANDLE CERTAIN PROBLEMS.
C      THE NEW VERSION OF THIS PACKAGE IS NOT INTENDED FOR USE ON
C      SUCH PROBLEMS AND ALL CODE RELATED TO IQUAD HAS THUS BEEN DELETED
C
C THE COMMON BLOCK /OPTION/ CONTAINS THE VARIABLES NOGAUS AND MAXDER.
C NOGAUS IS SET .EQ. 0 IN THE BLOCK DATA.  IT CAN BE CHANGED TO BE
C SET .EQ. 1 IF THE GAUSS-LEGENDRE COLLOCATION POINTS ARE NOT
C DESIRED WHEN NCC = 2 (SEE ABOVE AND COLPNT).  MAXDER IS SET
C .EQ. 5 IN THE BLOCK DATA AND ITS VALUE REPRESENTS THE
C MAXIMUM ORDER OF TIME INTEGRATION FORMULA ALLOWED.  ITS VALUE
C AFFECTS THE STORAGE REQUIRED IN WORK AND MAY BE CHANGED IF
C DESIRED.  SEE COSET FOR RESTRICTIONS.  THESE CHANGES MAY BE MADE BY
C THE USER BY ACCESSING /OPTION/ IN HIS CALLING PROGRAM (BEFORE THE
C FIRST CALL TO PDECOL) OR BY CHANGING THE DATA STATEMENT IN
C THE BLOCK DATA.
C
C-----------------------------------------------------------------------
C
C
C COMMUNICATION
C
C EACH SUBROUTINE IN THE PACKAGE CONTAINS A COMMUNICATION SUMMARY
C AS INDICATED BELOW.
C
C PACKAGE ROUTINES CALLED..   EVAL,INITAL,INTERP,STIFIB
C USER ROUTINES CALLED..      BNDRY
C CALLED BY..                 USERS MAIN PROGRAM
C FORTRAN FUNCTIONS USED..    ABS,MAX,FLOAT,SQRT
C-----------------------------------------------------------------------
      DIMENSION WORK(1),IWORK(1),XBKPT(1)
      INTEGER BOTPOS,SIZE
      COMMON /GEAR0/ DTUSED,NQUSED,NSTEP,NFE,NJE
      COMMON /GEAR1/ T,DTC,DTMN,DTMX,EPSC,UROUND,N,MFC,KFLAG,JSTART
      COMMON /GEAR9/ EPSJ,R0,IDUM(6)
      COMMON /OPTION/ NOGAUS,MAXDER
      COMMON /SIZES/ NIN,KOR,NC,NPD,NCPTS,NEQN,IFLAG
      COMMON/ABDPAR/NROWS,NCLBLK,NOVRLP,NELS,NELTOP,BOTPOS,SIZE,KRPT
      COMMON /ISTART/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,
     *                 IW11,IW12,IW13,IW14,IW15,IW16,IW17,IW18
      COMMON /IOUNIT/ LOUT
      IF (INDEX .EQ. 0) GO TO 60
      IF (INDEX .EQ. 2) GO TO 70
      IF (INDEX .EQ. 4) GO TO 80
      IF (INDEX .EQ. 3) GO TO 90
C-----------------------------------------------------------------------
C SEVERAL CHECKS ARE MADE HERE TO DETERMINE IF THE INPUT PARAMETERS
C HAVE LEGAL VALUES.  ERROR CHECKS ARE MADE ON INDEX, EPS, (T0-TOUT)*DT,
C DT, NINT, KORD, NCC, NPDE, MF, WHETHER THE BREAKPOINT SEQUENCE IS
C STRICTLY INCREASING, AND WHETHER THERE IS SUFFICIENT STORAGE
C PROVIDED FOR WORK AND IWORK.  PROBLEM DEPENDENT PARAMETERS ARE
C CALCULATED AND PLACED IN COMMON.
C-----------------------------------------------------------------------
      IERID = -6
      IF (INDEX .NE. 1) GO TO 320
      IERID = IERID - 1
      IF (EPS .LE. 0.) GO TO 320
      IERID = IERID - 1
      IF ((T0-TOUT)*DT .GT. 0.) GO TO 320
      IERID = IERID - 1
      IF (DT .EQ. 0.0) GO TO 320
      IERID = IERID - 1
      NIN = NINT
      IF (NIN .LT. 1) GO TO 320
      IERID = IERID - 1
      KOR = KORD
      IF (KOR .LT. 3  .OR.  KOR .GT. 20) GO TO 320
      IERID = IERID - 1
      NC = NCC
      IF (NCC .LT. 2  .OR.  NCC .GE. KOR) GO TO 320
      IERID = IERID - 1
      NPD = NPDE
      NPDE2 = NPD*NPD
      IF (NPDE .LT. 1) GO TO 320
      IERID = IERID - 1
      IF (MF.NE.22.AND.MF.NE.21.AND.MF.NE.12.AND.MF.NE.11) GO TO 320
      IERID = IERID - 1
      DO 10 K=1,NIN
        IF(XBKPT(K) .GE. XBKPT(K+1)) GO TO 320
   10 CONTINUE
      NCPTS = KOR + (NIN - 1) * (KOR - NCC)
      NEQN = NPDE * NCPTS
      IWSAVE = IWORK(1)
      IISAVE = IWORK(2)
      IW1 = 4*NPDE2 + 1
      IW2 = IW1 + 3*KORD*NCPTS
      IW3 = IW2 + NCPTS + KORD
      IW4 = IW3 + NCPTS
      IW5 = IW4 + NEQN
      IW6 = IW5 + NEQN
      IW7 = IW6 + NEQN
      IW8 = IW7 + NEQN
      IW9 = IW8 + NEQN
      IW10 = IW9 + 3*NPDE
      IW11 = IW10 + NEQN*(MAXDER+1)
      IW12 = IW11 + NPDE2
      IW13 = IW12 + NPDE2
      IW14 = IW13 + NPDE2
      IW15 = IW14 + NPDE2
      IW16 = IW15 + NPDE2
      IW17 = IW16 + NPDE
      IW18 = NCPTS + 1
      IERID = IERID - 1
      KRPT = KORD - NCC
      NROWS = KRPT*NPDE
      NELS = NROWS*NPDE*KORD
      NOVRLP = NCC*NPDE
      NCLBLK = KORD*NPDE
      NELTOP = NOVRLP*NPDE
      BOTPOS = NELTOP + NELS*NINT
      SIZE = BOTPOS + NELTOP
      IWSTOR = IW17 + SIZE - 1
      IISTOR = IW18 + NEQN - 1
      IF ( IWSAVE .LT. IWSTOR  .OR.  IISAVE .LT. IISTOR ) GO TO 335
C-----------------------------------------------------------------------
C PERFORM INITIALIZATION TASKS.
C-----------------------------------------------------------------------
      CALL INITAL(KOR,WORK(IW1),WORK(IW6),XBKPT,WORK(IW2),WORK(IW3),
     *            WORK(IW17),IWORK(IW18),IWORK)
      IF (IFLAG .NE. 0) GO TO 280
C-----------------------------------------------------------------------
C IF INITIAL VALUES OF CMAX OTHER THAN THOSE SET BELOW ARE DESIRED,
C THEY SHOULD BE SET HERE.  ALL CMAX(I) MUST BE POSITIVE.
C HAVING PROPER VALUES OF CMAX FOR THE PROBLEM BEING SOLVED IS AS
C IMPORTANT AS CHOOSING EPS (SEE ABOVE), SINCE ERRORS ARE
C MEASURED RELATIVE TO CMAX.  IF VALUES FOR DTMN OR DTMX, THE
C BOUNDS ON ABS(DT), OTHER THAN THOSE BELOW ARE DESIRED, THEY
C SHOULD BE SET BELOW.
C-----------------------------------------------------------------------
      DO 50 I = 1,NEQN
        I1 = I - 1
        WORK(IW4+I1) = ABS(WORK(IW6+I1))
        IF (WORK(IW4+I1) .LT. 1.) WORK(IW4+I1) = 1.
   50   WORK(IW10+I1) = WORK(IW6+I1)
      N = NEQN
      T = T0
      DTC = DT
      DTMN = ABS(DT)
      DTUSED = 0.
      EPSC = EPS
      MFC = MF
      JSTART = 0
      EPSJ =  SQRT(UROUND)
      NHCUT = 0
      KFLAG = 0
      TOUTP = T0
      IF ( T0 .EQ. TOUT ) GO TO 360
   60 DTMX = ABS(TOUT-TOUTP)*10.
      GO TO 140
C
   70 DTMX = ABS(TOUT-TOUTP)*10.
      IF ((T-TOUT)*DTC .GE. 0.) GO TO 340
      GO TO 150
C
   80 IF ((T-TOUT)*DTC .GE. 0.) GO TO 300
      JSTART = -1
      EPSC = EPS
      MFC = MF
      GO TO 100
C
   90 DTMX = DT
  100 IF ((T+DTC) .EQ. T) WRITE(LOUT,110)
  110 FORMAT(36H WARNING..  T + DT = T ON NEXT STEP.)
C-----------------------------------------------------------------------
C TAKE A TIME STEP BY CALLING THE INTEGRATOR.
C-----------------------------------------------------------------------
      CALL STIFIB (NEQN,WORK(IW10),WORK(IW4),WORK(IW5),WORK(IW6),
     *       WORK(IW7),WORK(IW8),WORK(IW17),IWORK(IW18),WORK,IWORK)
C
      KGO = 1 - KFLAG
      GO TO (120, 160, 220, 260, 280), KGO
C KFLAG  =   0,  -1,  -2,  -3   -4
C
  120 CONTINUE
C-----------------------------------------------------------------------
C NORMAL RETURN FROM INTEGRATOR.
C
C THE WEIGHTS CMAX(I) ARE UPDATED.  IF DIFFERENT VALUES ARE DESIRED,
C THEY SHOULD BE SET HERE.  A TEST IS MADE FOR EPS BEING TOO SMALL
C FOR THE MACHINE PRECISION.
C
C ANY OTHER TESTS OR CALCULATIONS THAT ARE REQUIRED AFTER EVERY
C STEP SHOULD BE INSERTED HERE.
C
C IF INDEX = 3, SAVE1 IS SET TO THE CURRENT C VALUES ON RETURN.
C IF INDEX = 2, DT IS CONTROLLED TO HIT TOUT (WITHIN ROUNDOFF
C ERROR), AND THEN THE CURRENT C VALUES ARE PUT IN SAVE1 ON RETURN.
C FOR ANY OTHER VALUE OF INDEX, CONTROL RETURNS TO THE INTEGRATOR
C UNLESS TOUT HAS BEEN REACHED.  THEN INTERPOLATED VALUES OF C ARE
C COMPUTED AND STORED IN SAVE1 ON RETURN.
C IF INTERPOLATION IS NOT DESIRED, THE CALL TO INTERP SHOULD BE
C REMOVED AND CONTROL TRANSFERRED TO STATEMENT 340 INSTEAD OF 360.
C-----------------------------------------------------------------------
      D = 0.
      DO 130 I = 1,NEQN
        I1 = I - 1
        AYI = ABS(WORK(IW10+I1))
        WORK(IW4+I1) = MAX(WORK(IW4+I1), AYI)
  130   D = D + (AYI/WORK(IW4+I1))**2
      D = D*(UROUND/EPS)**2
      IF (D .GT.  FLOAT(NEQN)) GO TO 240
      IF (INDEX .EQ. 3) GO TO 340
      IF (INDEX .EQ. 2) GO TO 150
  140 IF ((T-TOUT)*DTC .LT. 0.) GO TO 100
      CALL INTERP(TOUT,WORK(IW10),NEQN,WORK(IW6))
      GO TO 360
C
  150 IF (((T+DTC)-TOUT)*DTC .LE. 0.) GO TO 100
      IF (ABS(T-TOUT) .LE. 100.*UROUND*DTMX) GO TO 340
      IF ((T-TOUT)*DTC .GE. 0.) GO TO 340
      DTC = (TOUT - T)*(1. - 4.*UROUND)
      JSTART = -1
      GO TO 100
C-----------------------------------------------------------------------
C ON AN ERROR RETURN FROM INTEGRATOR, AN IMMEDIATE RETURN OCCURS IF
C KFLAG = -2 OR -4, AND RECOVERY ATTEMPTS ARE MADE OTHERWISE.
C TO RECOVER, DT AND DTMN ARE REDUCED BY A FACTOR OF .1 UP TO 10
C TIMES BEFORE GIVING UP.
C-----------------------------------------------------------------------
  160 WRITE (LOUT,170) T
  170 FORMAT(//35H KFLAG = -1 FROM INTEGRATOR AT T = ,E16.8/
     *   40H  ERROR TEST FAILED WITH ABS(DT) = DTMIN/)
  180 IF (NHCUT .EQ. 10) GO TO 200
      NHCUT = NHCUT + 1
      DTMN = .1*DTMN
      DTC = .1*DTC
      WRITE (LOUT,190) DTC
  190 FORMAT(25H  DT HAS BEEN REDUCED TO ,E16.8,
     *   26H  AND STEP WILL BE RETRIED//)
      JSTART = -1
      GO TO 100
C
  200 WRITE (LOUT,210)
  210 FORMAT(//44H PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT//)
      GO TO 340
C
  220 WRITE (LOUT,230) T,DTC
  230 FORMAT(//35H KFLAG = -2 FROM INTEGRATOR AT T = ,E16.8,6H  DT =,
     *  E16.8/52H  THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED//)
      GO TO 340
C
  240 WRITE (LOUT,250) T
  250 FORMAT(//37H INTEGRATION HALTED BY DRIVER AT T = ,E16.8/
     *   56H  EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION/)
      KFLAG = -2
      GO TO 340
C
  260 WRITE (LOUT,270) T
  270 FORMAT(//35H KFLAG = -3 FROM INTEGRATOR AT T = ,E16.8/
     *   45H  CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED/)
      GO TO 180
C
  280 WRITE (LOUT,290)
  290 FORMAT(//28H SINGULAR MATRIX ENCOUNTERED,
     *         35H PROBABLY DUE TO STORAGE OVERWRITES//)
      KFLAG = -4
      GO TO 340
C
  300 WRITE(LOUT,310) T,TOUT,DTC
  310 FORMAT(//45H INDEX = -1 ON INPUT WITH (T-TOUT)*DT .GE. 0./
     *   4H T =,E16.8,9H   TOUT =,E16.8,8H   DTC =,E16.8/
     *   44H INTERPOLATION WAS DONE AS ON NORMAL RETURN./
     *   41H DESIRED PARAMETER CHANGES WERE NOT MADE.)
      CALL INTERP(TOUT,WORK(IW10),NEQN,WORK(IW6))
      INDEX = -5
      RETURN
C
  320 WRITE(LOUT,330) IERID
  330 FORMAT(//24H ILLEGAL INPUT...INDEX= ,I3//)
      INDEX = IERID
      RETURN
C
  335 WRITE(LOUT,336) IWSTOR,IWSAVE,IISTOR,IISAVE
  336 FORMAT(//21H INSUFFICIENT STORAGE/24H  WORK MUST BE OF LENGTH,
     * I10,5X,12HYOU PROVIDED,I10/24H IWORK MUST BE OF LENGTH,I10,5X,
     * 12HYOU PROVIDED,I10//)
      INDEX = IERID
      RETURN
C
  340 TOUT = T
      DO 350 I = 1,NEQN
        I1 = I - 1
  350   WORK(IW6+I1) = WORK(IW10+I1)
  360 INDEX = KFLAG
      TOUTP = TOUT
      DT = DTUSED
      IF (KFLAG .NE. 0) DT = DTC
      RETURN
      END
      SUBROUTINE ADDA(PW,A,BC,NPDE)
*
      IMPLICIT REAL*8(A-H,O-Z)
*
C-----------------------------------------------------------------------
C SUBROUTINE ADDA ADDS ELEMENTS OF THE MATRIX A(Y,T) TO THE MATRIX
C STORED IN PW. THE ELEMENTS OF A(Y,T) ARE STORED IN THE INPUT PARAMETER
C MATRIX A(K,I,J). THIS ROUTINE EXTRACTS THE INFORMATION FROM THE LATTER
C ADDS IT TO PW IN AN APPROPRIATE FASHION, CREATING THE EFFECT OF ADDING
C A(Y,T) TO PW. ALSO, BOUNDARY CONDITION INFORMATION IS ADDED TO  THE FI
C AND LAST BLOCKS OF PW.
C
C DOCUMENTATION OF PARAMETERS
C
C PW   - THE JACOBIAN MATRIX TO WHICH THE MATRIX A(Y,T) IS ADDED.
C        IT IS A VECTOR OF DIMENSION SIZE. PW CONTAINS, COLUMN-WISE,
C        THE BLOCK, TOP, ASSOCIATED WITH THE LEFT BOUNDARY CONDITIONS,
C        FOLLOWED BY THE NINT BLOCKS ASSOCIATED WITH COLLOCATION POINTS
C        ON EACH SUBINTERVAL, FOLLOWED BY THE BLOCK, BOT, ASSOCIATED WIT
C        THE RIGHT BOUNDARY CONDITIONS. THE TOP BLOCK CONSISTS OF A SING
C        ROW OF TWO SUBBLOCKS, EACH OF SIZE NPDE BY NPDE. THE BOT BLOCK
C        THE SAME STRUCTURE. EACH OF THE OTHER BLOCKS IS A KRPT BY KORD
C        OF SUBBLOCKS, EACH OF SIZE NPDE BY NPDE. THE (I,J)-TH SUCH SUBB
C        CORRESPONDS TO THE J-TH NON-ZERO BASIS FUNCTION OF THE  SUBINTE
C        EVALUATED AT THE I-TH COLLOCATION POINT OF THAT SUBINTERVAL.
C A    - BASIS FUNCTION VALUES AT COLLOCATION POINTS. A IS USUALLY
C        A THREE DIMENSIONAL ARRAY OF SIZE (KORD,3,NCPTS), BUT IS
C        VIEWED HERE AS A ONE DIMENSIONAL ARRAY OF SIZE KORD*NCPTS*3.
C        A(K,J,I) CONTAINS THE VALUE OF THE (J-1)ST DERIVATIVE (J=1,2,3)
C        OF THE K-TH NON-ZERO BASIS FUNCTION (K=1,...,KORD) AT THE I-TH
C        COLLOCATION POINT.
C BC   - BOUNDARY CONDITION INFORMATION . BC IS USUALLY A THREE DIMENSIO
C        MATRIX OF SIZE (NPDE,NPDE,4) BUT HERE IT IS VIEWED AS A VECTOR
C        OF LENGTH 4*NPDE**2.
C NPDE - THE NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS.
C
C PACKAGE ROUTINES CALLED..  NONE
C USER ROUTINES CALLED..     NONE
C CALLED BY..                DIFFUN,PSETIB
C FORTRAN FUNCTIONS USED..   NONE
C
*
      DIMENSION PW(1),A(1),BC(NPDE,NPDE,4)
*
      INTEGER BOTPOS,SIZE
*
      COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPTS,NEQN,IDUM
      COMMON/ABDPAR/NROWS,NCLBLK,NOVRLP,NELS,NELTOP,BOTPOS,SIZE,KRPT
*
C     DOCUMENTATION OF VARIABLES IN COMMON :
C
C   NINT  - NUMBER OF SUBINTERVALS.
C   KORD  - NUMBER OF NON-ZERO BASIS FUNCTIONS PER SUBINTERVAL.
C   NCC   - NUMBER OF CONTINUITY CONDITIONS AT EACH INTERNAL MESH POINT.
C   NPDE  - NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS.
C   NCPTS - NUMBER OF COLLOCATION POINTS.
C   NEQN  - NUMBER OF DISCRETE EQUATIONS, NPDE * NCPTS
C
C   NROWS  - THE NUMBER OF ROWS IN EACH BLOCK TO BE PROVIDED TO CRDCMP.
C   NCLBLK - THE NUMBER OF COLUMNS IN EACH BLOCK TO BE PROVIDED TO CRDCM
C   NOVRLP - THE NUMBER OF COLUMNS OF OVERLAP BETWEEN THE BLOCKS PROVIDE
C            TO CRDCMP.
C   NELS   - THE NUMBER OF ELEMENTS IN ONE COLLOCATION BLOCK OF PW.
C   NELTOP - THE NUMBER OF ELEMENTS IN THE TOP BLOCK STORED IN PW.
C   BOTPOS - THE POSITION IN PW WHERE STORAGE LOCATIONS FOR THE INFORMAT
C            FOR THE BOTTOM BLOCK BEGIN.
C   SIZE   - THE LENGTH OF THE VECTOR PW.
C   KRPT   - THE NUMBER OF COLLOCATION POINTS PER SUBINTERVAL.
C-----------------------------------------------------------------------
C
C ADD THE BOUNDARY CONDITION INFORMATION, PROVIDED IN BC, TO THE TOP
C AND BOTTOM BLOCK ROWS OF PW.
C
      NPDE2 = NPDE*NPDE
*
      DO 200 J = 1,NPDE
*
        J1TOP = (J-1)*NPDE
*
C       J1TOP+1 IS THE POSITION IN PW OF THE FIRST ROW OF THE FIRST NPDE
C       SUBBLOCK OF THE TOP BLOCK.
*
        J2TOP = (J-1)*NPDE + NPDE2
*
C       J2TOP+1 IS THE POSITION IN PW OF THE FIRST ROW OF THE SECOND NPD
C       NPDE SUBBLOCK OF THE TOP BLOCK.
*
        J1BOT = BOTPOS + J1TOP
*
C       J1BOT+1 IS THE POSITION IN PW OF THE FIRST ROW OF THE FIRST NPDE
C       SUBBLOCK OF THE BOT BLOCK.
*
        J2BOT = BOTPOS + J2TOP
*
C       J2BOT+1 IS THE POSITION IN PW OF THE FIRST ROW OF THE SECOND NPD
C       NPDE SUBBLOCK OF THE BOT BLOCK.
*
        DO 100 I = 1,NPDE
*
          PW(J1TOP+I) = PW(J1TOP+I) + BC(I,J,1)
          PW(J2TOP+I) = PW(J2TOP+I) + BC(I,J,2)
          PW(J1BOT+I) = PW(J1BOT+I) + BC(I,J,3)
          PW(J2BOT+I) = PW(J2BOT+I) + BC(I,J,4)
*
  100   CONTINUE
*
  200 CONTINUE
*
C UPDATE THE REMAINING ROWS OF PW BY ADDING THE APPROPRIATE VALUES
C IN A TO PW.
*
      DO 600 I = 1,NINT
*
          IPOS = (I-1)*NELS + NELTOP
*
C         IPOS+1 IS THE POSITION IN PW OF THE START OF THE I-TH BLOCK
*
          IC = (I-1)*KRPT
*
C         IC+1 IS THE POSITION IN XC OF THE FIRST COLLOCATION POINT OF T
C         I-TH SUBINTERVAL
*
              DO 500 IP = 1,KRPT
*
                IPPOS = IPOS + (IP-1)*NPDE
*
C               IPPOS+1 IS THE POSITION IN PW CORRESONDING TO THE ROW OF
C               SUBBLOCKS ASSOCIATED WITH THE IP-TH COLLOCATION POINT OF
C               I-TH SUBINTERVAL
*
                IK = (IC+IP)*KORD*3
*
C               IK IS THE POSITION IN A(K,I,J) OF THE BASIS FUNCTION
C               INFORMATION CORRESPONDING TO THE IP-TH COLLOCATION POINT
C               THE I-TH SUBINTERVAL
*
                    DO 400 J = 1,KORD
*
                      JK = IK + J
*
C                     JK IS THE POSITION IN A(K,I,J) OF THE BASIS FUNCTI
C                     INFORMATION CORRESPONDING TO THE J-TH NON-ZERO BAS
C                     FUNCTION ASSOCIATED WITH THE IP-TH COLLOCATION POI
C                     OF THE I-TH SUBINTERVAL
*
                      JPOS = IPPOS + (J-1)*NROWS*NPDE
*
C                     JPOS+1 IS THE POSITION IN PW OF THE J-TH SUBBLOCK
C                     ROW OF SUBBLOCKS CORREPONDING TO THE IP-TH COLLOCA
C                     POINT OF THE I-TH SUBINTERVAL
*
                            DO 300 JJ = 1,NPDE
*
                              JJPOS = JPOS + (JJ-1)*NROWS + JJ
*
C                             JJPOS IS THE POSITION IN PW OF THE JJ-TH D
C                             ENTRY OF THE J-TH SUBBLOCK ABOVE.
*
C                             HERE WE ASSIGN THE VALUE OF THE J-TH NON-Z
C                             BASIS FUNCTION EVALUATED AT THE IP-TH COLL
C                             POINT OF THE I-TH SUBINTERVAL TO THE DIAGO
C                             OF THE (IP,J)-TH SUBBLOCK OF THE COLLOCATI
C                             BLOCK ASSOCIATED WITH THE I-TH SUBINTERVAL
*
                              PW(JJPOS) = PW(JJPOS) + A(JK)
*
  300                       CONTINUE
*
  400               CONTINUE
*
  500         CONTINUE
*
  600 CONTINUE
*
      RETURN
      END
      SUBROUTINE DIFFUN (N, T, Y, YDOT, IFLAG, PW, IPIV, WORK, IWORK)
*
      IMPLICIT REAL*8(A-H,O-Z)
*
C-----------------------------------------------------------------------
C STIFIB FINDS THE SOLUTION TO A SYSTEM OF O.D.E.'S OF THE FORM
C A(Y,T) (DY/DT) = G(Y,T). THIS REQUIRES THE CALCULATION OF THE QUANTITY
C A(Y,T)**-1 * G(Y,T). DIFFUN COMPUTES YDOT = A(Y,T)**-1 * G(Y,T) BY
C SETTING UP AND SOLVING THE SYSTEM A(Y,T) * YDOT = G(Y,T) FOR YDOT.
C
C DOCUMENTATION OF PARAMETERS
C
C N    - THE DIMENSION OF Y. (N = NEQN).
C T    - CURRENT VALUE OF THE INDEPENDENT VARIABLE T.
C Y    - THE VECTOR OF UNKNOWN BASIS FUNCTION COEFFICIENTS.
C        IT CONTAINS THE CURRENT APPROXIMATIONS OF THESE COEFFICIENTS.
C        THIS VECTOR IS ALSO KNOWN AS C IN OTHER PARTS OF THE PACKAGE.
C YDOT - TO BE COMPUTED AS INDICATED ABOVE.
C IFLAG- ERROR INDICATOR FROM THE MATRIX DECOMPOSITION ROUTINE.
C        IFLAG <> 0 IMPLIES FAILURE.
C PW   - COEFFICIENT MATRIX FOR DETERMINATION OF YDOT.
C IPIV - PIVOTING INFORMATION FROM THE MATRIX DECOMPOSITION ROUTINE.
C WORK - WORKING STORAGE VECTOR CONTAINING MANY OF THE IMPORTANT ARRAYS
C        USED IN THIS PACKAGE. COMPLETE DOCUMENTATION OF WORK IS GIVEN
C        IN THE PDECOL ROUTINE. HERE WE MENTION ONLY THOSE COMPONENTS
C        OF WORK THAT ARE USED IN THIS ROUTINE. EACH COMPONENT WILL BE
C        DESIGNATED BY ITS STARTING LOCATION IN WORK AND ALSO IT USUAL
C        NAME THROUGHOUT THE PACKAGE.
C
C        WORK(IW1), A      - BASIS FUNCTION VALUES AT COLLOCATION POINTS
C                            A IS USUALLY A THREE DIMENSIONAL ARRAY OF S
C                            (KORD,3,NCPTS), BUT IS VIEWED HERE AS A ONE
C                            DIMENSIONAL ARRAY OF SIZE KORD*NCPTS*3. A(K
C                            CONTAINS THE VALUE OF THE (J-1)ST DERIVATIV
C                            (J=1,2,3) OF THE K-TH NON-ZERO BASIS FUNCTI
C                            (K=1,...,KORD) AT THE I-TH COLLOCATION POIN
C        WORK(IW3), XC     - COLLOCATION POINTS.THIS VECTOR IS OF LENGTH
C        WORK(IW9), UVAL   - STORAGE FOR VALUES OF THE CURRENT SOLUTION
C                            APPROXIMATION, U(X), AND ITS FIRST AND SECO
C                            DERIVATIVES AT A SINGLE POINT, X.
C                            THIS IS A TWO DIMENSIONAL ARRAY OF SIZE (NP
C        WORK(IW14), DBDU  - BOUNDARY CONDITION INFORMATION. THIS TWO
C                            DIMENSIONAL ARRAY IS OF SIZE (NPDE,NPDE) AN
C                            CONTAINS THE FIRST DERIVATIVE WITH RESPECT
C                            OF THE LEFT HAND SIDE OF THE BOUNDARY CONDI
C                            EQUATIONS.
C
C        WORK(IW15), DBDUX - BOUNDARY CONDITION INFORMATION. THIS TWO
C                            DIMENSIONAL ARRAY IS OF SIZE (NPDE,NPDE) AN
C                            CONTAINS THE FIRST DERIVATIVE WITH RESPECT
C                            UX OF THE LEFT HAND SIDE OF THE BOUNDARY
C                            CONDITION EQUATIONS.
C        WORK(IW16), DZDT  - BOUNDARY CONDITION INFORMATION. THIS VECTOR
C                            LENGTH NPDE AND CONTAINS THE DERIVATIVE WIT
C                            RESPECT TO T OF THE RIGHT HAND SIDE OF THE
C                            BOUNDARY CONDITION EQUATIONS.
C
C IWORK - INTEGER WORKING STORAGE VECTOR. IT IS DOCUMENTED IN PDECOL.
C         HERE WE MENTION ONLY THE VECTOR IT CONTAINS THAT IS RELEVANT
C         TO THIS ROUTINE.
C
C         IWORK(1), ILEFT  - VECTOR OF BREAKPOINT SEQUENCE INFORMATION.
C                            THIS VECTOR IS OF LENGTH NCPTS.
C
C PACKAGE ROUTINES CALLED..  ADDA,DECB,GFUN,SOLB
C USER ROUTINES CALLED..     NONE
C CALLED BY..                STIFIB
C FORTRAN FUNCTIONS USED..   NONE
*
      DIMENSION Y(N),YDOT(N),PW(1),IPIV(1),WORK(1),IWORK(1)
*
      INTEGER BOTPOS,SIZE,NELTP1,BOTPS1
*
      COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IDUM0
      COMMON/ABDPAR/NROWS,NCLBLK,NOVRLP,NELS,NELTOP,BOTPOS,SIZE,KRPT
      COMMON /ISTART/ IW1,IDUM1,IW3,IDUM2(5),IW9,IDUM3(4),IW14,
     *                IW15,IW16,IDUM4(2)
*
C DOCUMENTATION OF VARIABLES IN COMMON :
C
C NINT  - NUMBER OF SUBINTERVALS.
C KORD  - NUMBER OF NON-ZERO BASIS FUNCTIONS PER SUBINTERVAL.
C NCC   - NUMBER OF CONTINUITY CONDITIONS AT EACH INTERNAL MESH POINT.
C NPDE  - NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS.
C NCPTS - NUMBER OF COLLOCATION POINTS.
C NEQN  - NUMBER OF DISCRETE EQUATIONS, NPDE * NCPTS
C
C NROWS  - THE NUMBER OF ROWS IN EACH BLOCK TO BE PROVIDED TO CRDCMP.
C NCLBLK - THE NUMBER OF COLUMNS IN EACH BLOCK TO BE PROVIDED TO CRDCMP.
C NOVRLP - THE NUMBER OF COLUMNS OF OVERLAP BETWEEN THE BLOCKS PROVIDED
C          TO CRDCMP.
C NELS   - THE NUMBER OF ELEMENTS IN ONE COLLOCATION BLOCK OF PW.
C NELTOP - THE NUMBER OF ELEMENTS IN THE TOP BLOCK STORED IN PW.
C BOTPOS - THE POSITION IN PW WHERE THE STORAGE LOCATIONS FOR THE INFORM
C          FOR THE BOTTOM BLOCK BEGIN.
C SIZE   - THE LENGTH OF THE VECTOR PW.
C KRPT   - THE NUMBER OF COLLOCATION POINTS PER SUBINTERVAL.
C
C THE COMMON BLOCK /ISTART/ CONTAINS THE STARTING LOCATIONS, IN WORK
C AND IWORK OF THE COMPONENT VECTORS CONTAINING THE ARRAYS OF
C INFORMATION USED THROUGHTOUT THIS PACKAGE.
*
C-----------------------------------------------------------------------
*
*
C EVALUATE THE FUNCTION G(Y,T). IT IS RETURNED IN THE VECTOR YDOT,
C AND IS USED AS THE RIGHT HAND SIDE OF THE NEWTON SYSTEM. BOUNDARY
C CONDITION INFORMATION IS ALSO UPDATED BY THIS CALL TO GFUN.
*
      CALL GFUN (T, Y, YDOT, NPDE, NCPTS, WORK(IW1), WORK, WORK(IW14),
     *           WORK(IW15), WORK(IW16), WORK(IW3), WORK(IW9), IWORK)
*
C THE FOLLOWING CALL TO GFUN IS EQUVIVALENT TO THE ABOVE CALL WITH THE
C ACTUAL PARAMETER NAMES SUBSTITUED FOR THE WORK VECTOR ADDRESSES.
C
C     CALL GFUN (T, Y, YDOT, NPDE, NCPTS, A, BC, DBDU,
C    *           DBDUX, DZDT, XC, UVAL, ILEFT)
*
C INITIALIZE PW.
*
      DO 10 I = 1,SIZE
*
        PW(I) = 0.0D+0
*
10    CONTINUE
*
C ADD INFORMATION FROM THE MATRIX A(K,I,J) TO THE MATRIX PW.
C PW WILL THEN CONTAIN THE MATRIX A(Y,T).
*
      CALL ADDA (PW, WORK(IW1), WORK, NPDE)
*
C THE FOLLOWING CALL TO ADDA IS EQUVIVALENT TO THE ABOVE CALL WITH THE
C ACTUAL PARAMETER NAMES SUBSTITUED FOR THE WORK VECTOR ADDRESSES.
C
C     CALL ADDA (PW, A, BC, NPDE)
*
C FACTOR THE MATRIX PW.
*
      BOTPS1 = BOTPOS + 1
      NELTP1 = NELTOP + 1
      IFLAG = 0
      print *,' CRDCMP, NEQN = ',neqn
      CALL CRDCMP(NEQN,PW,NPDE,NOVRLP,PW(NELTP1),NROWS,NCLBLK,NINT,
     *            PW(BOTPS1),NPDE,IPIV,IFLAG)
*
      IF (IFLAG .NE. 0) GO TO 1021
*
C SOLVE THE TRIANGULAR SYSTEMS RETURNING THE SOLUTION IN YDOT.
*
      print *,' CRSLVE, nEQN = ',neqn
      CALL CRSLVE(PW,NPDE,NOVRLP,PW(NELTP1),NROWS,NCLBLK,NINT,
     *            PW(BOTPS1),NPDE,IPIV,YDOT,YDOT)
*
 1021 CONTINUE
      RETURN
      END
      SUBROUTINE INITAL(K,A,RHS,X,XT,XC,PW,IPIV,ILEFT)
*
      IMPLICIT REAL*8(A-H,O-Z)
*
C-----------------------------------------------------------------------
C INITAL IS CALLED ONLY ONCE BY PDECOL TO PERFORM INITIALIZATION TASKS.
C THESE TASKS INCLUDE - 1) DEFINING THE PIECEWISE POLYNOMIAL SPACE
C BREAKPOINT SEQUENCE, 2) CALLING THE SUBROUTINE COLPNT TO DEFINE THE
C REQUIRED COLLOCATION POINTS, 3) DEFINING THE PIECEWISE POLYNOMIAL SPAC
C BASIS FUNCTION VALUES (PLUS FIRST AND SECOND DERIVATIVE VALUES) AT
C THE COLLOCATION POINTS, AND 4) DEFINING THE INITIAL BASIS FUNCTION
C COEFFICIENTS WHICH DETERMINE THE PIECEWISE POLYNOMIAL WHICH
C INTERPOLATES THE USER SUPPLIED INITIAL CONDITIONS AT THE
C COLLOCATION POINTS. THIS LAST STEP INVOLVES THE CONSTRUCTION OF A LINE
C SYSTEM OF THE FORM PW*C = RHS, WHERE C IS THE VECTOR OF
C LENGTH NPDE*NCPTS CONTAINING THE UNKNOWN COEFFICIENTS TO BE DETERMINED
C BY THIS ROUTINE. RHS CONSISTS OF NCPTS VECTORS EACH OF LENGTH NPDE.
C THE I-TH SUCH VECTOR IS THE RIGHT HAND SIDE OF THE INITIAL CONDITIONS
C EVALUATED AT THE I-TH COLLOCATION POINT. THE MATRIX PW CONSISTS OF AN
C NCPTS BY NCPTS TWO DIMENSIONAL ARRAY OF BLOCKS, WITH EACH BLOCK BEING
C SIZE NPDE BY NPDE. THE (I,J)-TH SUCH BLOCK CONSISTS OF THE IDENTITY MA
C MULTIPLIED BY THE J-TH BASIS FUNCTION EVALUATED AT THE I-TH COLLOCATIO
C POINT. BECAUSE OF SPECIAL PROPERTIES OF THE BASIS FUNCTIONS, ONLY KORD
C OF THEM ARE NON-ZERO WHEN EVALUATED AT THE SET OF COLLOCATION POINTS
C ASSOCIATED WITH A SINGLE SUBINTERVAL. THIS PROPERTY RESULTS IN THE PW
C HAVING AN ALMOST BLOCK DIAGONAL STRUCTURE.
C
C DOCUMENTATION OF PARAMETERS
C
C K     - EQUAL TO KORD, THE NUMBER OF NON-ZERO BASIS FUNCTIONS
C         PER SUBINTERVAL; USED FOR DIMENSIONING PURPOSES.
C A     - BASIS FUNCTION VALUES AT COLLOCATION POINTS. A IS
C         A THREE DIMENSIONAL ARRAY OF SIZE (KORD,3,NCPTS).
C         A(K,J,I) CONTAINS THE VALUE OF THE (J-1)ST DERIVATIVE (J=1,2,3
C         OF THE K-TH NON-ZERO BASIS FUNCTION (K=1,...,KORD) AT THE I-TH
C         COLLOCATION POINT.
C RHS   - IN THIS ROUTINE THE RIGHT HAND SIDE OF THE INITAL CONDITIONS
C         EVALUATED AT EACH OF THE NCPTS COLLOCATION POINTS IS ASSIGNED
C         TO RHS.  ON OUTPUT THIS VECTOR CONTAINS THE COEFFICIENTS WHICH
C         DETERMINE THE PIECE-WISE POLYNOMIAL INTERPOLATING THE INITIAL
C         CONDITIONS AT THE COLLOCATION POINTS.  IT IS OF SIZE NPDE*NCPT
C X     - THE MESHPOINT SEQUENCE. IT IS OF SIZE NINT + 1.
C XT    - THE BREAKPOINT SEQUENCE. IT IS OF SIZE NCPTS + KORD.
C XC    - THE COLLOCATION POINT SEQUENCE. IT IS OF SIZE NCPTS.
C PW    - IN INITAL, THIS VECTOR, WHICH IS OF DIMENSION SIZE, WILL BE
C         ASSIGNED VALUES AS INDICATED ABOVE. PW CONTAINS, COLUMN-WISE,
C         THE BLOCK, TOP, ASSOCIATED WITH THE LEFT BOUNDARY CONDITIONS,
C         FOLLOWED BY THE NINT BLOCKS ASSOCIATED WITH COLLOCATION POINTS
C         ON EACH SUBINTERVAL, FOLLOWED BY THE BLOCK, BOT, ASSOCIATED WI
C         THE RIGHT BOUNDARY CONDITIONS. THE TOP BLOCK CONSISTS OF A SIN
C         ROW OF TWO SUBBLOCKS, EACH OF SIZE NPDE BY NPDE. THE BOT BLOCK
C         THE SAME STRUCTURE. EACH OF THE OTHER BLOCKS IS A KRPT BY KORD
C         OF SUBBLOCKS, EACH OF SIZE NPDE BY NPDE. THE (I,J)-TH SUCH SUB
C         CORRESPONDS TO THE J-TH NON-ZERO BASIS FUNCTION OF THE  SUBINT
C         EVALUATED AT THE I-TH COLLOCATION POINT OF THAT SUBINTERVAL.
C         ON OUTPUT PW CONTAINS THE FACTORED COEFFICIENT MATRIX
C         USED TO DETERMINE THE COMPONENTS OF RHS.
C IPIV  - PIVOTING INFORMATION FROM THE FACTORIZATION OF PW.
C         IT IS OF SIZE NPDE*NCPTS.
C ILEFT - BREAKPOINT INFORMATION. IT IS OF SIZE NCPTS.
C
C PACKAGE ROUTINES CALLED..  BSPLVD,COLPNT,CRDCMP,INTERV,CRSLVE
C USER ROUTINES CALLED..     UINIT
C CALLED BY..                PDECOL
C FORTRAN FUNCTIONS USED..   MAX0,MIN0
*
      DIMENSION A(K,3,1),RHS(1),X(1),XT(1),XC(1),PW(1),IPIV(1),ILEFT(1)
      INTEGER BOTPOS,SIZE,BOTPS1
*
*
      COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IFLAG
      COMMON/ABDPAR/NROWS,NCLBLK,NOVRLP,NELS,NELTOP,BOTPOS,SIZE,KRPT
*
C     DOCUMENTATION OF VARIABLES IN COMMON :
C
C   NINT  - NUMBER OF SUBINTERVALS.
C   KORD  - NUMBER OF NON-ZERO BASIS FUNCTIONS PER SUBINTERVAL.
C   NCC   - NUMBER OF CONTINUITY CONDITIONS AT EACH INTERNAL MESH POINT.
C   NPDE  - NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS.
C   NCPTS - NUMBER OF COLLOCATION POINTS.
C   NEQN  - NUMBER OF DISCRETE EQUATIONS, NPDE * NCPTS
C   IFLAG - ERROR INDICATOR FOR LINEAR SYSTEM FACTORIZATION ROUTINE.
C
C   NROWS  - THE NUMBER OF ROWS IN EACH BLOCK TO BE PROVIDED TO CRDCMP.
C   NCLBLK - THE NUMBER OF COLUMNS IN EACH BLOCK TO BE PROVIDED TO CRDCM
C   NOVRLP - THE NUMBER OF COLUMNS OF OVERLAP BETWEEN THE BLOCKS PROVIDE
C            TO CRDCMP.
C   NELS   - THE NUMBER OF ELEMENTS IN ONE COLLOCATION BLOCK OF PW.
C   NELTOP - THE NUMBER OF ELEMENTS IN THE TOP BLOCK STORED IN PW.
C   BOTPOS - THE POSITION IN PW WHERE STORAGE LOCATIONS FOR THE INFORMAT
C            FOR THE BOTTOM BLOCK BEGIN.
C   SIZE   - THE LENGTH OF THE VECTOR PW.
C   KRPT   - THE NUMBER OF COLLOCATION POINTS PER SUBINTERVAL.
C
C-----------------------------------------------------------------------
*
C SET UP THE PIECEWISE POLYNOMIAL SPACE BREAKPOINT SEQUENCE.
C FIRST SET UP THE SEQUENCE AT EACH ENDPOINT.
      DO 10 I=1,KORD
*
            XT(NCPTS+I) = X(NINT+1)
            XT(I) = X(1)
*
 10   CONTINUE
*
C SET UP THE REMAINING POINTS OF THE SEQUENCE.
*
      DO 30 I=2,NINT
*
        IK = (I-2)*KRPT + KORD
*
C       IK+1 IS THE POSITION IN XT OF THE START OF THE BREAKPOINT INFORM
C       ASSOCIATED WITH THE KRPT COLLOCATION POINTS LOCATED ON THE I-TH
C       SUBINTERVAL.
*
        DO 20 J=1,KRPT
*
            XT(IK+J) = X(I)
*
 20     CONTINUE
*
 30   CONTINUE
*
C SET UP COLLOCATION POINTS ARRAY XC.
*
       CALL COLPNT(X, XC, XT)
*
C GENERATE THE ILEFT ARRAY.  STORE THE BASIS FUNCTION VALUES IN THE
C ARRAY A.  THE ARRAY A IS DIMENSIONED A(KORD,3,NCPTS) AND A(K,J,I)
C CONTAINS THE VALUE OF THE (J-1)-ST DERIVATIVE (J = 1,2,3) OF THE K-TH
C NONZERO BASIS FUNCTION (K = 1, ... ,KORD) AT THE I-TH COLLOCATION
C POINT (I = 1, ... ,NCPTS).  SET UP RHS FOR INTERPOLATING THE INITIAL
C CONDITIONS AT THE COLLOCATION POINTS.  ASSIGN VALUES FROM THE ABOVE MA
C A, TO THE MATRIX PW.
C
C INITIALIZE PW TO ZERO.
C
      DO 1130 I = 1,SIZE
*
                PW(I) = 0.0D+0
*
 1130 CONTINUE
*
C FIRST, SET UP THE TOP BLOCK OF THE MATRIX, MAKING USE OF THE FACT
C THAT ONLY THE FIRST B-SPLINE HAS A NON-ZERO VALUE AT THE LEFT END
C POINT. ALSO, SET UP THE FIRST NPDE COMPONENTS OF THE RIGHT HAND
C SIDE.
*
      MFLAG = -2
       CALL INTERV(XT,NCPTS,XC(1),ILEFT(1),MFLAG)
*
C     HERE INTERV COMPUTES THE VALUE OF ILEFT(1). MFLAG IS INITIALIZED T
C     -2 TO INSTRUCT INTERV TO BEGIN ITS SEARCH FOR THE APPROPRIATE ILEF
C     WITH THE VALUE 1.
*
       CALL BSPLVD(XT,KORD,XC(1),ILEFT(1),A(1,1,1),3)
C
C     BSPLVD IS CALLED TO COMPUTE THE COMPONENTS OF A(K,I,J) ASSOCIATED
C     THE FIRST COLLOCATION POINT.
C
       CALL UINIT(XC(1),RHS(1),NPDE)
C
C     UINIT IS CALLED TO EVALUATE THE FIRST NPDE COMPONENTS OF THE RHS V
C     THE RIGHT HAND SIDE OF THE INITIAL CONDITIONS, EVALUATED AT THE FI
C     COLLOCATION POINT IS ASSIGNED TO THE RHS VECTOR.
*
      DO 4040 N = 1,NPDE
*
          NN = (N-1)*NPDE + N
          PW(NN) = A(1,1,1)
*
C         HERE WE ASSIGN A(1,1,1), THE FIRST BASIS FUNCTION EVALUATED AT
C         FIRST COLLOCATION POINT, TO THE DIAGONAL OF THE FIRST (I.E TOP
C         BLOCK OF PW. THE INDEX NN MOVES THROUGH THE LOCATIONS IN PW
C         CORRESPONDING TO THIS DIAGONAL.
*
 4040 CONTINUE
*
C THE NINT BLOCKS OF THE MAIN MATRIX WILL NOW BE SET UP. EACH SUCH BLOCK
C CONSISTS OF A KRPT BY KORD TWO DIMENSIONAL ARRAY OF SUBBLOCKS, EACH OF
C SIZE NPDE BY NPDE.  THE CORRESPONDING COMPONENTS OF THE RIGHT HAND SID
C VECTOR ARE ALSO ASSIGNED.
*
*
        DO 1040 I = 1,NINT
*
        IPOS = (I-1)*NELS + NELTOP
*
C       IPOS+1 IS THE POSITION IN PW OF THE START OF THE I-TH BLOCK
*
        IC = (I-1)*KRPT
*
C       IC+1 IS THE POSITION IN XC OF THE FIRST COLLOCATION POINT OF THE
C       I-TH SUBINTERVAL
*
C       FOR EACH COLLOCATION POINT OF THE I-TH SUBINTERVAL,
C       CONSTRUCT THE CORRESPONDING SUBBLOCK ROW OF THE I-TH BLOCK.
*
          DO 1030 IP = 1,KRPT
*
          IPPOS = IPOS + (IP-1)*NPDE
*
C         IPPOS+1 IS THE POSITION IN PW CORRESONDING TO THE ROW OF SUBBL
C         ASSOCIATED WITH THE IP-TH COLLOCATION POINT OF THE I-TH SUBINT
*
          ICC = IC + IP + 1
*
C         ICC IS THE ABSOLUTE POSITION IN XC OF THE IP-TH COLLKCATION PO
C         OF THE I-TH SUBINTERVAL
*
          IN = (ICC-1)*NPDE
*
C         IN+1 IS THE POSITION IN THE RHS VECTOR WHERE THE VALUES FOR TH
C         RIGHT HAND SIDE OF THE INITAL CONDITIONS, EVALUATED AT THE ICC
C         COLLOCATION POINT ARE STORED.
*
C         COMPUTE INFORMATION FOR ICC-TH COLLOCATION POINT
*
           CALL INTERV(XT,NCPTS,XC(ICC),ILEFT(ICC),MFLAG)
           CALL BSPLVD(XT,KORD,XC(ICC),ILEFT(ICC),A(1,1,ICC),3)
           CALL UINIT(XC(ICC),RHS(IN+1),NPDE)
*
C         THE ABOVE THREE CALLS PERFORM THE SAME TASKS AS MENTIONED ABOV
C         EXCEPT THAT ALL EVAULATIONS ARE DONE AT THE ICC-TH COLLOCATION
C         POINT INSTEAD OF AT THE FIRST.
*
            DO 1010 J = 1,KORD
*
C           FOR J RANGING OVER THE KORD SUBBLOCKS OF THE IP-TH SUBBLOCK
C           OF THE I-TH SUBINTERVAL CONSTRUCT THE CORRESPONDING COMPONEN
C           PW.
*
            JPOS = IPPOS + (J-1)*NROWS*NPDE
*
C           JPOS+1 IS THE POSITION IN PW OF THE J-TH SUBBLOCK OF THE ROW
C           SUBBLOCKS CORREPONDING TO THE IP-TH COLLOCATION POINT OF THE
C           SUBINTERVAL
*
                DO 1005 L = 1,NPDE
*
                LPOS = JPOS + (L-1)*NROWS + L
*
C               LPOS IS THE POSITION IN PW OF THE (L,L)-TH ENTRY OF THE
C               J-TH SUBBLOCK ABOVE.
*
                PW(LPOS) = A(J,1,ICC)
*
C               HERE WE ASSIGN THE J-TH NON-ZERO BASIS FUNCTION OF THE I
C               SUBINTERVAL, EVALUATED AT THE ICC-TH COLLOCATION POINT,
C               DIAGONAL OF THE (IP,J)-TH NPDE BY NPDE SUBBLOCK OF THE I
C               COLLOCATION BLOCK.
*
 1005           CONTINUE
*
 1010       CONTINUE
*
 1030    CONTINUE
*
 1040  CONTINUE
*
C NOW, SET UP THE BOTTOM BLOCK, USING THE FACT THAT ONLY THE LAST
C B-SPLINE BASIS FUNCTION IS NON-ZERO AT THE RIGHT END POINT.
C SIMULTANEOUSLY, SET UP THE CORRESPONDING PART OF THE RIGHT HAND SIDE.
*
       CALL INTERV(XT,NCPTS,XC(NCPTS),ILEFT(NCPTS),MFLAG)
       CALL BSPLVD(XT,KORD,XC(NCPTS),ILEFT(NCPTS),A(1,1,NCPTS),3)
      IN = NEQN - NPDE +1
       CALL UINIT(XC(NCPTS),RHS(IN),NPDE)
*
      DO 3040 N = 1,NPDE
*
        NN = ((N-1) + (NCC-1)*NPDE)*NPDE + N + BOTPOS
        PW(NN) = A(KORD,1,NCPTS)
*
 3040 CONTINUE
*
C LU DECOMPOSE THE MATRIX PW.
*
      NELTP1 = NELTOP + 1
      BOTPS1 = BOTPOS + 1
      IFLAG = 0
*
      print *,' CRDCMP, NEQN = ',neqn
       CALL CRDCMP(NEQN,PW,NPDE,NOVRLP,PW(NELTP1),NROWS,NCLBLK,NINT,
     *            PW(BOTPS1),NPDE,IPIV,IFLAG)
*
*
      IF (IFLAG.NE.0) GO TO 1021
*
C SOLVE THE LINEAR SYSTEM   PW*Z = RHS.  THIS GIVES THE BASIS FUNCTION
C COEFFICIENTS FOR THE INITIAL CONDITIONS. THE SOLUTION, Z, IS RETURNED
C IN THE VECTOR RHS.
*
       print *,' CRSLVE, NEQN = ',neqn
       CALL CRSLVE(PW,NPDE,NOVRLP,PW(NELTP1),NROWS,NCLBLK,NINT,
     *            PW(BOTPS1),NPDE,IPIV,RHS,RHS)
*
*
 1021 CONTINUE
*
      RETURN
      END
      SUBROUTINE PSETIB (C, PW, N0, CON, MITER, IFLAG, A, ILEFT, XC,
     *   UVAL, SAVE2,IPIV,CMAX,DFDU,DFDUX,DFDUXX,DZDT,DBDU,DBDUX,BC,
     *   NPDE)
*
      IMPLICIT REAL*8(A-H,O-Z)
*
C-----------------------------------------------------------------------
C STIFIB IS USED TO SOLVE A SYSTEM OF ODE'S OF THE FORM A(C,T)*DC/DT = G
C WHEN NEWTON'S METHOD IS APPLIED TO SOLVE THE DISCRETE NON-LINEAR SYSTE
C ARISING FROM THE ABOVE SYSTEM OF ODE'S, THE JACOBIAN MUST BE COMPUTED.
C THE MATRIX THAT MUST BE COMPUTED HAS THE FORM A(C,T) + CON*(DG/DC),
C WHERE THE CONSTANT, CON, IS DEFINED BELOW. PSETIB IS CALLED BY STIFIB
C CONSTRUCT AND FACTOR THIS MATRIX, WHICH IS TO BE STORED IN PW.
C DG/DC IS COMPUTED ACCORDING TO THE VALUE OF MITER DEFINED BELOW.
C FINALLY, PW IS SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER
C SOLUTION OF LINEAR SYSTEMS WITH PW AS COEFFICIENT MATRIX.
C
C DOCUMENTATION OF PARAMETERS
C
C C     - THE VECTOR OF UNKNOWN COEFFICIENTS; THE VECTOR C
C         CONTAINS THE CURRENT APPROXIMATE VALUES OF THESE COEFFICIENTS.
C         THIS VECTOR IS OF SIZE NCPTS*NPDE.
C
C PW    - IN PSETIB, THIS VECTOR, WHICH IS OF DIMENSION SIZE, WILL BE
C         ASSIGNED VALUES AS INDICATED ABOVE. PW CONTAINS, COLUMN-WISE,
C         THE BLOCK, TOP, ASSOCIATED WITH THE LEFT BOUNDARY CONDITIONS,
C         FOLLOWED BY THE NINT BLOCKS ASSOCIATED WITH COLLOCATION POINTS
C         ON EACH SUBINTERVAL, FOLLOWED BY THE BLOCK, BOT, ASSOCIATED WI
C         THE RIGHT BOUNDARY CONDITIONS. THE TOP BLOCK CONSISTS OF A SIN
C         ROW OF TWO SUBBLOCKS, EACH OF SIZE NPDE BY NPDE. THE BOT BLOCK
C         THE SAME STRUCTURE. EACH OF THE OTHER BLOCKS IS A KRPT BY KORD
C         OF SUBBLOCKS, EACH OF SIZE NPDE BY NPDE. THE (I,J)-TH SUCH SUB
C         CORRESPONDS TO THE J-TH NON-ZERO BASIS FUNCTION OF THE  SUBINT
C         EVALUATED AT THE I-TH COLLOCATION POINT OF THAT SUBINTERVAL.
C         ON OUTPUT PW CONTAINS THE FACTORED COEFFICIENT MATRIX.
C
C NO    - EQUAL TO NEQN. USED TO PROVIDE RUN-TIME
C         DIMENSIONING FOR SEVERAL ARRAYS.
C
C CON   - THIS CONSTANT CONTAINS THE PRODUCT OF TWO VALUES MENTIONED
C         ABOVE IN REFERENCE TO THE DEFINITION OF PW. IT EQUALS
C         -H * EL(1), WHERE EL(1) IS A COEFFICIENT FROM THE NUMERICAL
C         FORMULA USED TO DISCRETIZE THE SYSTEM OF ODE'S.
C
C MITER - THIS FLAG INDICATES WHETHER OR NOT AN ANALYTIC JACOBIAN IS
C         TO BE PROVIDED BY THE USER; MITER = 1 IMPLIES AN ANALYTIC
C         JACOBIAN WILL BE PROVIDED THROUGH THE ROUTINE DERIVF;
C         MITER = 2 IMPLIES THAT A DIVIDED-DIFFERENCE APPROXIMATION
C         TO THE JACOBIAN MUST BE COMPUTED.
C
C IFLAG - ERROR INDICATOR FROM THE MATRIX DECOMPOSITION ROUTINE.
C         IFLAG <> 0 IMPLIES A FAILURE.
C
C A     - BASIS FUNCTION VALUES AT COLLOCATION POINTS. A IS USUALLY
C         A THREE DIMENSIONAL ARRAY OF SIZE (KORD,3,NCPTS), BUT IS
C         VIEWED HERE AS A ONE DIMENSIONAL ARRAY OF SIZE KORD*NCPTS*3.
C         A(K,J,I) CONTAINS THE VALUE OF THE (J-1)ST DERIVATIVE (J=1,2,3
C         OF THE K-TH NON-ZERO BASIS FUNCTION (K=1,...,KORD) AT THE I-TH
C         COLLOCATION POINT.
C
C ILEFT - POINTERS TO BREAKPOINT SEQUENCE. THIS VECTOR IS OF SIZE NCPTS.
C
C XC    - COLLOCATION POINTS. THIS VECTOR IS OF LENGTH NCPTS.
C
C UVAL  - WORKING STORAGE FOR VALUES OF U, UX, AND UXX AT ONE POINT.
C         THIS IS A TWO DIMENSIONAL ARRAY OF SIZE (NPDE,3).
C
C SAVE2 - WORKING STORAGE FOR THE TIME INTEGRATION METHOD.
C         THIS VECTOR IS OF LENGTH NPDE*NCPTS. IT CONTAINS NCPTS VECTORS
C         EACH OF LENGTH NPDE. EACH OF THESE VECTORS CORRESPONDS TO AN
C         EVALUATION OF THE RIGHT HAND SIDE OF THE SYSTEM OF PDE'S AT
C         A COLLOCATION POINT.
C
C IPIV  - PIVOT INFORMATION FOR THE LU DECOMPOSED JACOBIAN MATRIX PW.
C         THIS VECTOR IS OF LENGTH NPDE*NCPTS.
C
C CMAX  - VALUES USED IN ESTIMATING TIME INTEGRATION ERRORS.
C         THIS INFORMATION IS USED BY THE ROUTINE DIFFF WHEN A NUMERICAL
C         APPROXIMATION TO THE JACOBIAN MUST BE COMPUTED.
C
C DFDU  - WORKING STORAGE USED TO COMPUTE THE JACOBIAN MATRIX.
C         THIS TWO DIMENSIONAL ARRAY IS OF SIZE (NPDE,NPDE) AND CONTAINS
C         THE DERIVATIVE WITH RESPECT TO U OF THE RIGHT HAND SIDE OF THE
C         SYSTEM OF PDE'S, EVALUATED AT A COLLOCATION POINT.
C
C DFDUX - WORKING STORAGE USED TO COMPUTE THE JACOBIAN MATRIX.
C         THIS TWO DIMENSIONAL ARRAY IS OF SIZE (NPDE,NPDE) AND CONTAINS
C         THE DERIVATIVE WITH RESPECT TO UX OF THE RIGHT HAND SIDE OF TH
C         SYSTEM OF PDE'S, EVALUATED AT A COLLOCATION POINT.
C
C DFDUXX- WORKING STORAGE USED TO COMPUTE THE JACOBIAN MATRIX.
C         THIS TWO DIMENSIONAL ARRAY IS OF SIZE (NPDE,NPDE) AND CONTAINS
C         THE DERIVATIVE WITH RESPECT TO UXX OF THE RIGHT HAND SIDE OF T
C         SYSTEM OF PDE'S, EVALUATED AT A COLLOCATION POINT.
C
C DZDT  - BOUNDARY CONDITION INFORMATION. THIS VECTOR IS OF LENGTH NPDE
C         AND CONTAINS THE DEPIVATIVE WITH RESPECT TO T OF THE RIGHT HAN
C         SIDE OF THE BOUNDARY CONDITION EQUATIONS.
C
C DBDU  - BOUNDARY CONDITION INFORMATION. THIS TWO DIMENSIONAL ARRAY IS
C         SIZE (NPDE,NPDE) AND CONTAINS THE FIRST DERIVATIVE WITH RESPEC
C         OF THE LEFT HAND SIDE OF THE BOUNDARY CONDITION EQUATIONS.
C
C DBDUX - BOUNDARY CONDITION INFORMATION. THIS TWO DIMENSIONAL ARRAY IS
C         SIZE (NPDE,NPDE) AND CONTAINS THE FIRST DERIVATIVE WITH RESPEC
C         OF THE LEFT HAND SIDE OF THE BOUNDARY CONDITION EQUATIONS.
C
C BC    - BOUNDARY CONDITION INFORMATION. BC IS USUALLY A THREE DIMENSIO
C         MATRIX OF SIZE (NPDE,NPDE,4) BUT HERE IT IS VIEWED AS A VECTOR
C         LENGTH 4*NPDE**2.
C
C
C NPDE - NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS
C
C PACKAGE ROUTINES CALLED..  ADDA,CRDCMP,DIFFF,EVAL,GFUN
C USER ROUTINES CALLED..     BNDRY,DERIVF
C CALLED BY..                STIFIB
C FORTRAN FUNCTIONS USED..   ABS,FLOAT,MAX,MIN,SQRT
C
      INTEGER BOTPOS,SIZE,NELTP1,BOTPS1
      DIMENSION PW(SIZE),C(1),CMAX(1)
      DIMENSION A(1),ILEFT(1),BC(1),XC(1),UVAL(NPDE,3),SAVE2(1),IPIV(1)
      DIMENSION DFDU(NPDE,NPDE),DFDUX(NPDE,NPDE),DFDUXX(NPDE,NPDE)
      DIMENSION DZDT(NPDE),DBDU(NPDE,NPDE),DBDUX(NPDE,NPDE)
*
      COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPTS,NEQN,IDUM
      COMMON /GEAR1/ T,H,DUM1(3),UROUND,N,IDUM1(3)
      COMMON /GEAR9/ DUM2,R0,IDUM2(6)
      COMMON/ABDPAR/NROWS,NCLBLK,NOVRLP,NELS,NELTOP,BOTPOS,SIZE,KRPT
*
C DOCUMENTATION OF VARIABLES IN COMMON :
C
C NINT  - NUMBER OF SUBINTERVALS.
C KORD  - NUMBER OF NON-ZERO BASIS FUNCTIONS PER SUBINTERVAL.
C NCC   - NUMBER OF CONTINUITY CONDITIONS AT EACH INTERNAL MESH POINT.
C NPDE  - NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS.
C NCPTS - NUMBER OF COLLOCATION POINTS.
C NEQN  - NUMBER OF DISCRETE EQUATIONS, NPDE * NCPTS
C
C T      - THE INDEPENDENT VARIABLE REPRESENTING THE TIME DIMENSION.
C H      - THE CURRENT STEPSIZE BEING USED IN THE STIFIB ROUTINE.
C UROUND - THE UNIT ROUND-OFF OF THE MACHINE.
C N      - THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS
C          IN THE SYSTEMS TO BE SOLVED BY STIFIB.
C
C R0     - SCALING FACTOR - USED BY THE DIFFF ROUTINE AND CALCULATED IN
C          PSETIB.
C
C NROWS  - THE NUMBER OF ROWS IN EACH BLOCK TO BE PROVIDED TO CRDCMP.
C NCLBLK - THE NUMBER OF COLUMNS IN EACH BLOCK TO BE PROVIDED TO CRDCMP.
C NOVRLP - THE NUMBER OF COLUMNS OF OVERLAP BETWEEN THE BLOCKS PROVIDED
C          TO CRDCMP.
C NELS   - THE NUMBER OF ELEMENTS IN ONE COLLOCATION BLOCK OF PW.
C NELTOP - THE NUMBER OF ELEMENTS IN THE TOP BLOCK STORED IN PW.
C BOTPOS - THE POSITION IN PW WHERE THE INFORMATION FOR THE BOTTOM
C          BLOCK IS STORED.
C SIZE   - THE LENGTH OF THE VECTOR PW.
C KRPT   - THE NUMBER OF COLLOCATION POINTS PER SUBINTERVAL.
C
C
C-----------------------------------------------------------------------
*
C     INITIALIZE PW.
*
      DO 10 I=1,SIZE
*
        PW(I) = 0.0D+0
*
 10   CONTINUE
*
C     IF THE USER DOES NOT PROVIDE THE ANALYTIC JACOBIAN, (I.E.
C     MITER IS NOT EQUAL TO 1), THEN, BY CALLING GFUN, COMPUTE
C     THE REQUIRED INFORMATION TO BE USED LATER, DURING THE CALL
C     TO DIFFF, IN THE CALCULATION OF AN APPROXIMATION TO THE JACOBIAN.
*
      IF ( MITER .EQ. 1 ) GO TO 25
*
      CALL GFUN (T, C, SAVE2, NPDE, NCPTS,A,BC,DBDU,DBDUX,DZDT,XC,
     *           UVAL,ILEFT)
*
C     COMPUTE R0 FOR USE IN THE DIFFF ROUTINE.
*
      D = 0.0D+0
*
      DO 20 I = 1,N
*
        D = D + SAVE2(I)**2
*
 20   CONTINUE
*
      R0 =  ABS(H)* SQRT(D/FLOAT(N0))*1.0D+3*UROUND
*
   25 CONTINUE
*
C     FOR EACH OF THE NINT SUBINTERVALS COMPUTE THE CORRESPONDING
C     BLOCK OF THE COEFFICIENT MATRIX PW.
*
      DO 300 I=1,NINT
*
        IPOS = (I-1)*NELS + NELTOP
*
C       IPOS+1 IS THE POSITION IN PW OF THE START OF THE I-TH BLOCK
*
        IC = (I-1)*KRPT
*
C       IC+1 IS THE POSITION IN XC OF THE FIRST COLLOCATION POINT OF THE
C       I-TH SUBINTERVAL
*
C       FOR EACH COLLOCATION POINT OF THE I-TH SUBINTERVAL,
C       CONSTRUCT THE CORRESPONDING SUBBLOCK ROW OF THE I-TH BLOCK.
*
        DO 250 IP = 1,KRPT
*
          IPPOS = IPOS + (IP-1)*NPDE
*
C         IPPOS+1 IS THE POSITION IN PW CORRESONDING TO THE START OF THE
C         OF SUBBLOCKS ASSOCIATED WITH THE IP-TH COLLOCATION POINT OF TH
C         SUBINTERVAL
*
          IK = (IC+IP)*KORD*3
*
C         IK IS THE POSITION IN A OF THE BASIS FUNCTIGN INFORMATION CORR
C         ING TO THE IP-TH COLLOCATION POINT OF THE I-TH SUBINTERVAL
*
          ICC = IC + IP + 1
*
C         ICC IS THE ABSOLUTE POSITION IN XC OF THE IP-TH COLLKCATION PO
C         OF THE I-TH SUBINTERVAL
*
          CALL EVAL(ICC,NPDE,C,UVAL,A,ILEFT)
*
C         EVAL RETURNS AN EVALUATION OF THE SOLUTION COMPONENTS AND THEI
C         FIRST AND SECOND DERIVATIVES IN THE MATRIX UVAL, WITH THE EVAL
C         UATIONS BEING DONE AT THE ICC-TH COLLOCATION POINT.
*
          IF ( MITER .EQ. 1 )
     *       CALL DERIVF(T,XC(ICC),UVAL,UVAL(1,2),UVAL(1,3),
     *                              DFDU,DFDUX,DFDUXX,NPDE)
*
C         DERIVF RETURNS AN EVALUATION OF THE JACOBIAN OF F AS WELL
C         AS THE DERIVATIVES WITH RESPECT TO THE FIRST AND SECOND
C         DERIVATIVES OF THE SOLUTION IN THE MATRICES DFDU, DFDUX,
C         AND DFDUXX. THE EVALUATION IS DONE AT THE ICC-TH COLLOCATION
C         POINT.
*
          IF ( MITER .EQ. 2 )
     *       CALL DIFFF(T,XC(ICC),ICC,UVAL,UVAL(1,2),UVAL(1,3),
     *                     DFDU,DFDUX,DFDUXX,NPDE,CMAX,SAVE2)
*
C         IN THIS CASE THE JACOBIAN IS CALCULATED NUMERICALLY BY THE DIF
C         ROUTINE
*
          DO 200 J = 1,KORD
*
C           FOR J RANGING OVER THE KORD SUBBLOCKS OF THE IP-TH SUBBLOCK
C           OF THE I-TH SUBINTERVAL CONSTRUCT THE CORRESPONDING COMPONEN
C           PW. (THESE CORRESPOND TO THE KORD NON-ZERO BASIS FUNCTIONS O
C           I-TH SUBINTERVAL).
*
            JK = IK + J
*
C           JK IS THE POSITION IN A OF THE BASIS FUNCTION INFORMATION
C           CORRESPONDING TO THE J-TH NON-ZERO BASIS FUNCTION ASSOCIATED
C           WITH THE IP-TH COLLOCATION POINT OF THE I-TH SUBINTERVAL
*
            JK2 = JK + KORD
*
C           JK2 IS THE POSITION IN A OF THE FIRST DERIVATIVE INFORMATION
C           FOR THE ABOVE BASIS FUNCTION
*
            JK3 = JK2 + KORD
*
C           JK3 IS THE POSITION IN A OF THE SECOND DERIVATIVE INFORMATIO
C           FOR THE ABOVE BASIS FUNCTION
*
            JPOS = IPPOS + (J-1)*NROWS*NPDE
*
C           JPOS+1 IS THE POSITION IN PW OF THE START OF THE J-TH SUBBLO
C           THE ROW OF SUBBLOCKS CORREPONDING TO THE IP-TH COLLOCATION P
C           OF THE I-TH SUBINTERVAL
*
            DO 150 JJ = 1,NPDE
*
              DO 100 KK = 1,NPDE
*
                JJPOS = JPOS + (JJ-1)*NROWS + KK
*
C               JJPOS IS THE POSITION IN PW OF THE (JJ,KK)-TH ENTRY OF T
C               J-TH SUBBLOCK ABOVE.
*
C               WE NOW ASSIGN THE APPROPRIATE COMPONENT OF (DG/DC), IN T
C               OF DF/DU, DF/DUX, AND DF/DUXX AND A(K,I,J) TO PW.
*
                PW(JJPOS) = DFDU(KK,JJ)*A(JK) + DFDUX(KK,JJ)*A(JK2)
     *                        + DFDUXX(KK,JJ)*A(JK3)
*
  100         CONTINUE
*
  150       CONTINUE
*
  200     CONTINUE
*
  250   CONTINUE
*
  300 CONTINUE
*
C     THE NEW VERSION DOES NOT ALLOW COLLOCATION AT THE ENDPOINTS SINCE
C     THE SIZES OF THE TOP AND BOTTOM BLOCKS WOULD MAKE THEM UNSUITABLE
C     FOR USE IN COLROW. THUS AT THIS POINT THERE IS NO NEED TO MAKE
C     ADJUSTMENTS TO THE ENTRIES OF PW BASED ON WHETHER OR NOT SOME
C     COMPONENTS OF THE SYSTEM OF PDE'S HAVE NO BOUNDARY CONDITIONS IMPO
C     ON THEM, AS IS DONE IS THE ORIGINAL VERSION.
*
C     AT THIS POINT PW CONTAINS (DG/DC). THE REMAINING STEPS ARE TO MULT
C     BY THE CONSTANT, CON, AND THEN ADD ON THE MATRIX A(Y,T)
*
      DO 85 I = 1,SIZE
*
        PW(I) = PW(I)*CON
*
   85 CONTINUE
*
C ADD MATRIX A(C,T) TO PW.
*
        CALL ADDA (PW, A, BC, NPDE)
*
C DO LU DECOMPOSITION ON PW.
*
      NELTP1 = NELTOP + 1
      BOTPS1 = BOTPOS + 1
      IFLAG = 0
*
      print *,' CRDCMP, NEQN = ',neqn
      CALL CRDCMP (NEQN, PW, NPDE, NOVRLP, PW(NELTP1), NROWS, NCLBLK,
     *             NINT,PW(BOTPS1), NPDE, IPIV, IFLAG)
*
      RETURN
      END
      SUBROUTINE STIFIB (N0,Y,YMAX,ERROR,SAVE1,SAVE2,SAVE3,
     *                   PW,IPIV,WORK,IWORK)
      IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C STIFIB PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE
C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM,
C   A(Y,T)*(DY/DT) = G(Y,T),   WHERE   Y = (Y(1),Y(2), ... ,Y(N)).
C STIFIB IS FOR USE WHEN THE MATRICES A AND DG/DY HAVE AN ALMOST BLOCK
C DIAGONAL FORM.  THE DEPENDENCE OF A(Y,T) ON Y IS ASSUMED TO BE WEAK.
C
C---------------------------------------------------------------------
C THE MODIFICATION TO THE ORIGINAL STIFIB PROGRAM CONSISTS ENTIRELY OF T
C REMOVAL OF THE BAND LINEAR EQUATION SOLVER, SOLB, AND ITS SUBSTITUTION
C BY CRSLVE, THE ALMOST BLOCK DIAGONAL LINEAR EQUATION SOLVER. THE CHANG
C IN THE STIFIB ROUTINE ARE LIMITED TO:
C
C       (1) REPLACEMENT OF THE CALL TO SOLV BY A CALL TO CRSLVE
C       (2) THE REMOVAL OF SOME VARIABLES USED SOLELY TO PROVIDE
C           INFORMATION TO SOLB AND THE INSERTION OF NEW VARIABLES
C           TO PASS INFORMATION TO CRSLVE.
C-----------------------------------------------------------------------
C
C REFERENCE
C
C   HINDMARSH, A.C., PRELIMINARY DOCUMENTATION OF GEARIB.. SOLUTION
C     OF IMPLICIT SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS WITH
C     BANDED JACOBIANS, LAWRENCE LIVERMORE LAB, UCID-30130, FEBRUARY
C     1976.
C
C COMMUNICATION WITH STIFIB IS DONE WITH THE FOLLOWING VARIABLES..
C
C   Y       AN N0 BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES
C             AND THEIR SCALED DERIVATIVES.  LMAX IS 13 FOR THE ADAMS
C             METHODS AND 6 FOR THE GEAR METHODS.  LMAX - 1 = MAXDER
C             IS THE MAXIMUM ORDER AVAILABLE.  SEE SUBROUTINE COSET.
C             Y(I,J+1) CONTAINS THE J-TH DERIVATIVE OF Y(I), SCALED BY
C             H**J/FACTORIAL(J)  (J = 0,1,...,NQ).
C   N0      A CONSTANT INTEGER .GE. N, USED FOR DIMENSIONING PURPOSES.
C   T       THE INDEPENDENT VARIABLE. T IS UPDATED ON EACH STEP TAKEN.
C   H       THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP.
C             H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE
C             PROBLEM.  H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS
C             SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM.
C   HMIN,   THE MINIMUM AND MAXIMUM ABSOLUTE VALUE OF THE STEP SIZE
C    HMAX     TO BE USED FOR THE STEP.  THESE MAY BE CHANGED AT ANY
C             TIME, BUT WILL NOT TAKE EFFECT UNTIL THE NEXT H CHANGE.
C   EPS     THE RELATIVE ERROR BOUND.  SEE DESCRIPTION IN PDECOL.
C   UROUND  THE UNIT ROUNDOFF OF THE MACHINE.
C   N       THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS.
C   MF      THE METHOD FLAG.  SEE DESCRIPTION IN PDECOL.
C   KFLAG   A COMPLETION CODE WITH THE FOLLOWING MEANINGS..
C                     0  THE STEP WAS SUCCESFUL.
C                    -1  THE REQUESTED ERROR COULD NOT BE ACHIEVED
C                          WITH ABS(H) = HMIN.
C                    -2  THE REQUESTED ERROR IS SMALLER THAN CAN
C                          BE HANDLED FOR THIS PROBLEM.
C                    -3  CORRECTOR CONVERGENCE COULD NOT BE
C                          ACHIEVED FOR ABS(H) = HMIN.
C                    -4  SINGULAR A-MATRIX ENCOUNTERED.
C             ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF T AND
C             THE Y ARRAY ARE AS OF THE BEGINNING OF THE LAST
C             STEP, AND H IS THE LAST STEP SIZE ATTEMPTED.
C   JSTART  AN INTEGER USED ON INPUT AND OUTPUT.
C             ON INPUT, IT HAS THE FOLLOWING VALUES AND MEANINGS..
C                     0  PERFORM THE FIRST STEP.
C                 .GT.0  TAKE A NEW STEP CONTINUING FROM THE LAST.
C                 .LT.0  TAKE THE NEXT STEP WITH A NEW VALUE OF
C                          H, EPS, N, AND/OR MF.
C             ON EXIT, JSTART IS NQ, THE CURRENT ORDER OF THE METHOD.
C   YMAX    AN ARRAY OF N ELEMENTS WITH WHICH THE ESTIMATED LOCAL
C             ERRORS IN Y ARE COMPARED.
C   ERROR   AN ARRAY OF N ELEMENTS.  ERROR(I)/TQ(2) IS THE ESTIMATED
C             ONE-STEP ERROR IN Y(I).
C   SAVE1,SAVE2,SAVE3   THREE WORKING STORAGE ARRAYS, EACH OF LENGTH N.
C   PW      A BLOCK OF LOCATIONS USED FOR THE CHORD ITERATION
C             MATRIX.  SEE DESCRIPTION IN PDECOL.
C   IPIV    AN INTEGER ARRAY OF LENGTH N FOR PIVOT INFORMATION.
C   WORK,IWORK   WORKING ARRAYS WHICH ARE USED TO PASS APPROPRIATE
C                  ARRAYS TO OTHER SUBROUTINES. SEE DESCRIPTION IN
C                  PDECOL.
C
C PACKAGE ROUTINES CALLED..  COSET,DIFFUN,PSETIB,RES,CRSLVE
C USER ROUTINES CALLED..     NONE
C CALLED BY..                PDECOL
C FORTRAN FUNCTIONS USED..   ABS,MAX,MIN,FLOAT
C-----------------------------------------------------------------------
      DIMENSION Y(N0,1),YMAX(N0),ERROR(N0),SAVE1(N0),SAVE2(N0),
     *          SAVE3(N0),PW(1),IPIV(1),WORK(1),IWORK(1)
*
      COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD
      COMMON /ISTART/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,IW11,
     *                IW12,IW13,IW14,IW15,IW16,IW17,IW18
      COMMON /GEAR1/ T,H,HMIN,HMAX,EPS,UROUND,N,MF,KFLAG,JSTART
      COMMON /GEAR0/ HUSED,NQUSED,NSTEP,NFE,NJE
      COMMON /OPTION/ NOGAUS,MAXDER
      COMMON/ABDPAR/NROWS,NCLBLK,NOVRLP,NELS,NELTOP,BOTPOS,SIZE,KRPT
*
      INTEGER BOTPOS,SIZE,NELTP1,BOTPS1
      DIMENSION EL(13),TQ(4)
      DATA EL(2)/1./, OLDL0/1./, TQ(1)/0./, IER/0/
      KFLAG = 0
      TOLD = T
      IF (JSTART .GT. 0) GO TO 200
      IF (JSTART .NE. 0) GO TO 120
C-----------------------------------------------------------------------
C ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL YDOT IS
C CALCULATED.  RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED
C IN A SINGLE STEP.  IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL
C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10.  IF A FAILURE
C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2
C FOR THE NEXT INCREASE.
C-----------------------------------------------------------------------
      NQ = 1
      IER = 0
      CALL DIFFUN (N, T, Y, SAVE1, IER, PW, IPIV, WORK, IWORK)
      IF ( IER .NE. 0 ) GO TO 685
      DO 110 I = 1,N
  110   Y(I,2) = H*SAVE1(I)
      METH = MF/10
      MITER = MF - 10*METH
      L = 2
      IDOUB = 3
      RMAX = 1.D+4
      RC = 0.
      CRATE = 1.
      EPSOLD = EPS
      HOLD = H
      MFOLD = MF
      NOLD = N
      NSTEP = 0
      NSTEPJ = 0
      NFE = 0
      NJE = 1
      IRET = 3
      GO TO 130
C-----------------------------------------------------------------------
C IF THE CALLER HAS CHANGED METH, COSET IS CALLED TO SET
C THE COEFFICIENTS OF THE METHOD.  IF THE CALLER HAS CHANGED
C N, EPS, OR METH, THE CONSTANTS E, EDN, EUP, AND BND MUST BE RESET.
C E IS A COMPARISON FOR ERRORS OF THE CURRENT ORDER NQ. EUP IS
C TO TEST FOR INCREASING THE ORDER, EDN FOR DECREASING THE ORDER.
C BND IS USED TO TEST FOR CONVERGENCE OF THE CORRECTOR ITERATES.
C IF THE CALLER HAS CHANGED H, Y MUST BE RESCALED.
C IF H OR METH HAS BEEN CHANGED, IDOUB IS RESET TO L + 1 TO PREVENT
C FURTHER CHANGES IN H FOR THAT MANY STEPS.
C-----------------------------------------------------------------------
  120 IF (MF .EQ. MFOLD) GO TO 150
      MEO = METH
      MIO = MITER
      METH = MF/10
      MITER = MF - 10*METH
      MFOLD = MF
      IF (MITER .NE. MIO) IWEVAL = MITER
      IF (METH .EQ. MEO) GO TO 150
      IDOUB = L + 1
      IRET = 1
  130 CALL COSET (METH, NQ, EL, TQ)
      LMAX = MAXDER + 1
      RC = RC*EL(1)/OLDL0
      OLDL0 = EL(1)
  140 FN =  FLOAT(N)
      EDN = FN*(TQ(1)*EPS)**2
      E   = FN*(TQ(2)*EPS)**2
      EUP = FN*(TQ(3)*EPS)**2
      BND = FN*(TQ(4)*EPS)**2
      GO TO (160, 170, 200), IRET
  150 IF ((EPS .EQ. EPSOLD) .AND. (N .EQ. NOLD)) GO TO 160
      EPSOLD = EPS
      NOLD = N
      IRET = 1
      GO TO 140
  160 IF (H .EQ. HOLD) GO TO 200
      RH = H/HOLD
      H = HOLD
      IREDO = 3
      GO TO 175
  170 RH = MAX(RH,HMIN/ ABS(H))
  175 RH = MIN(RH,HMAX/ ABS(H),RMAX)
      R1 = 1.
      DO 180 J = 2,L
        R1 = R1*RH
        DO 180 I = 1,N
  180     Y(I,J) = Y(I,J)*R1
      H = H*RH
      RC = RC*RH
      IDOUB = L + 1
      IF (IREDO .EQ. 0) GO TO 690
C-----------------------------------------------------------------------
C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY
C MULTIPLYING THE Y ARRAY BY THE PASCAL TRIANGLE MATRIX.
C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT  H*EL(1).
C WHEN RC DIFFERS FROM 1 BY MORE THAN 30 PERCENT, OR THE CALLER HAS
C CHANGED MITER, IWEVAL IS SET TO MITER TO FORCE PW TO BE UPDATED.
C IN ANY CASE, PW IS UPDATED AT LEAST EVERY 40-TH STEP.
C PW IS THE CHORD ITERATION MATRIX A - H*EL(1)*(DG/DY).
C-----------------------------------------------------------------------
  200 IF ( ABS(RC-1.) .GT. 0.3) IWEVAL = MITER
      IF (NSTEP .GE. NSTEPJ+40) IWEVAL = MITER
      T = T + H
      DO 210 J1 = 1,NQ
        DO 210 J2 = J1,NQ
          J = (NQ + J1) - J2
          DO 210 I = 1,N
  210       Y(I,J) = Y(I,J) + Y(I,J+1)
C-----------------------------------------------------------------------
C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN.  A CONVERGENCE TEST IS
C MADE ON THE R.M.S. NORM OF EACH CORRECTION, USING BND, WHICH
C IS DEPENDENT ON EPS.  THE SUM OF THE CORRECTIONS IS ACCUMULATED
C IN THE VECTOR ERROR(I).  THE Y ARRAY IS NOT ALTERED IN THE CORRECTOR
C LOOP.  THE UPDATED Y VECTOR IS STORED TEMPORARILY IN SAVE1.
C THE UPDATED H*YDOT IS STORED IN SAVE2.
C-----------------------------------------------------------------------
  220 DO 230 I = 1,N
        SAVE2(I) = Y(I,2)
  230   ERROR(I) = 0.
      M = 0
      CALL RES (T, H, Y, SAVE2, SAVE3, NPDE, NCPTS, WORK(IW1), IWORK,
     *  WORK, WORK(IW14), WORK(IW15), WORK(IW16), WORK(IW3), WORK(IW9))
      NFE = NFE + 1
      IF (IWEVAL .LE. 0) GO TO 350
C-----------------------------------------------------------------------
C IF INDICATED, THE MATRIX PW IS REEVALUATED BEFORE STARTING THE
C CORRECTOR ITERATION.  IWEVAL IS SET TO 0 AS AN INDICATOR
C THAT THIS HAS BEEN DONE.  PW IS COMPUTED AND PROCESSED IN PSETIB.
C-----------------------------------------------------------------------
      IWEVAL = 0
      RC = 1.
      NJE = NJE + 1
      NSTEPJ = NSTEP
      CON = -H*EL(1)
      CALL PSETIB (Y, PW, N0, CON, MITER, IER, WORK(IW1), IWORK,
     *     WORK(IW3),WORK(IW9),SAVE2,IPIV,YMAX,WORK(IW11),WORK(IW12),
     *     WORK(IW13),WORK(IW16),WORK(IW14),WORK(IW15),WORK,NPDE)
      IF (IER .NE. 0) GO TO 420
C-----------------------------------------------------------------------
C COMPUTE THE CORRECTOR ERROR, R SUB M, AND SOLVE THE LINEAR SYSTEM
C WITH THAT AS RIGHT-HAND SIDE AND PW AS COEFFICIENT MATRIX,
C USING THE LU DECOMPOSITION OF PW.
C-----------------------------------------------------------------------
  350 CONTINUE
      NELTP1 = NELTOP + 1
      BOTPS1 = BOTPOS + 1
*
       print *,' CRSLVE, NEQN = ',neqn
      CALL CRSLVE( PW, NPDE, NOVRLP, PW(NELTP1), NROWS, NCLBLK,
     *             NINT, PW(BOTPS1), NPDE, IPIV, SAVE3, SAVE3)
*
  370 D = 0.
      DO 380 I = 1,N
        ERROR(I) = ERROR(I) + SAVE3(I)
        D = D + (SAVE3(I)/YMAX(I))**2
        SAVE1(I) = Y(I,1) + EL(1)*ERROR(I)
  380   SAVE2(I) = Y(I,2) + ERROR(I)
C-----------------------------------------------------------------------
C TEST FOR CONVERGENCE.  IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE
C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST.
C-----------------------------------------------------------------------
  400 IF (M .NE. 0) CRATE = MAX(.9*CRATE,D/D1)
      IF ((D*MIN(1.D-0,2.*CRATE)) .LE. BND) GO TO 450
      D1 = D
      M = M + 1
      IF (M .EQ. 3) GO TO 410
      CALL RES(T, H, SAVE1, SAVE2, SAVE3, NPDE, NCPTS, WORK(IW1), IWORK,
     *  WORK, WORK(IW14), WORK(IW15), WORK(IW16), WORK(IW3), WORK(IW9))
      GO TO 350
C-----------------------------------------------------------------------
C THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES.
C IF THE MATRIX PW IS NOT UP TO DATE, IT IS REEVALUATED FOR THE
C NEXT TRY.  OTHERWISE THE Y ARRAY IS RETRACTED TO ITS VALUES
C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE.  IF NOT, A
C NO-CONVERGENCE EXIT IS TAKEN.
C-----------------------------------------------------------------------
  410 NFE = NFE + 2
      IF (IWEVAL .EQ. -1) GO TO 440
  420 T = TOLD
      RMAX = 2.
      DO 430 J1 = 1,NQ
        DO 430 J2 = J1,NQ
          J = (NQ + J1) - J2
          DO 430 I = 1,N
  430       Y(I,J) = Y(I,J) - Y(I,J+1)
      IF ( ABS(H) .LE. HMIN*1.00001) GO TO 680
      RH = .25
      IREDO = 1
      GO TO 170
  440 IWEVAL = MITER
      GO TO 220
C-----------------------------------------------------------------------
C THE CORRECTOR HAS CONVERGED.  IWEVAL IS SET TO -1 TO SIGNAL
C THAT PW MAY NEED UPDATING ON SUBSEQUENT STEPS.  THE ERROR TEST
C IS MADE AND CONTROL PASSES TO STATEMENT 500 IF IT FAILS.
C-----------------------------------------------------------------------
 450  IWEVAL = -1
      NFE = NFE + M
      D = 0.
      DO 460 I = 1,N
  460   D = D + (ERROR(I)/YMAX(I))**2
      IF (D .GT. E) GO TO 500
C-----------------------------------------------------------------------
C AFTER A SUCCESSFUL STEP, UPDATE THE Y ARRAY.
C CONSIDER CHANGING H IF IDOUB = 1.  OTHERWISE DECREASE IDOUB BY 1.
C IF IDOUB IS THEN 1 AND NQ .LT. MAXDER, THEN ERROR IS SAVED FOR
C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP.
C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER
C BY ONE IS CONSIDERED ALSO.  A CHANGE IN H IS MADE ONLY IF IT IS BY A
C FACTOR OF AT LEAST 1.1.  IF NOT, IDOUB IS SET TO 10 TO PREVENT
C TESTING FOR THAT MANY STEPS.
C-----------------------------------------------------------------------
      KFLAG = 0
      IREDO = 0
      NSTEP = NSTEP + 1
      HUSED = H
      NQUSED = NQ
      DO 470 J = 1,L
        DO 470 I = 1,N
  470     Y(I,J) = Y(I,J) + EL(J)*ERROR(I)
      IF (IDOUB .EQ. 1) GO TO 520
      IDOUB = IDOUB - 1
      IF (IDOUB .GT. 1) GO TO 700
      IF (L .EQ. LMAX) GO TO 700
      DO 490 I = 1,N
  490   Y(I,LMAX) = ERROR(I)
      GO TO 700
C-----------------------------------------------------------------------
C THE ERROR TEST FAILED.  KFLAG KEEPS TRACK OF MULTIPLE FAILURES.
C RESTORE T AND THE Y ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE
C TO TRY THE STEP AGAIN.  COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR
C ONE LOWER ORDER.
C-----------------------------------------------------------------------
  500 KFLAG = KFLAG - 1
      T = TOLD
      DO 510 J1 = 1,NQ
        DO 510 J2 = J1,NQ
          J = (NQ + J1) - J2
          DO 510 I = 1,N
  510       Y(I,J) = Y(I,J) - Y(I,J+1)
      RMAX = 2.
      IF ( ABS(H) .LE. HMIN*1.00001) GO TO 660
      IF (KFLAG .LE. -3) GO TO 640
      IREDO = 2
      PR3 = 1.D+20
      GO TO 540
C-----------------------------------------------------------------------
C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS
C PR1, PR2, AND PR3 ARE COMPUTED, BY WHICH H COULD BE DIVIDED
C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY.
C IN THE CASE OF FAILURE, PR3 = 1.E20 TO AVOID AN ORDER INCREASE.
C THE SMALLEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN
C ACCORDINGLY.  IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE
C ADDITIONAL SCALED DERIVATIVE.
C-----------------------------------------------------------------------
  520 PR3 = 1.D+20
      IF (L .EQ. LMAX) GO TO 540
      D1 = 0.
      DO 530 I = 1,N
  530   D1 = D1 + ((ERROR(I) - Y(I,LMAX))/YMAX(I))**2
      ENQ3 = .5/ FLOAT(L+1)
      PR3 = ((D1/EUP)**ENQ3)*1.4 + 1.4D-6
  540 ENQ2 = .5/ FLOAT(L)
      PR2 = ((D/E)**ENQ2)*1.2 + 1.2D-6
      PR1 = 1.D+20
      IF (NQ .EQ. 1) GO TO 560
      D = 0.
      DO 550 I = 1,N
  550   D = D + (Y(I,L)/YMAX(I))**2
      ENQ1 = .5/ FLOAT(NQ)
      PR1 = ((D/EDN)**ENQ1)*1.3 + 1.3D-6
  560 IF (PR2 .LE. PR3) GO TO 570
      IF (PR3 .LT. PR1) GO TO 590
      GO TO 580
  570 IF (PR2 .GT. PR1) GO TO 580
      NEWQ = NQ
      RH = 1./PR2
      GO TO 620
  580 NEWQ = NQ - 1
      RH = 1./PR1
      GO TO 620
  590 NEWQ = L
      RH = 1./PR3
      IF (RH .LT. 1.1) GO TO 610
      DO 600 I = 1,N
  600   Y(I,NEWQ+1) = ERROR(I)*EL(L)/ FLOAT(L)
      GO TO 630
  610 IDOUB = 10
      GO TO 700
  620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1)) GO TO 610
C-----------------------------------------------------------------------
C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS.
C IN ANY CASE H IS RESET ACCORDING TO RH AND THE Y ARRAY IS RESCALED.
C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE.
C-----------------------------------------------------------------------
      IF (NEWQ .EQ. NQ) GO TO 170
  630 NQ = NEWQ
      L = NQ + 1
      IRET = 2
      GO TO 130
C-----------------------------------------------------------------------
C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURED.
C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE
C Y ARRAY HAVE ERRORS OF THE WRONG ORDER.  HENCE THE FIRST
C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1.  THEN
C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED.
C AFTER A TOTAL OF 7 FAILURES, AN EXIT IS TAKEN WITH KFLAG = -2.
C-----------------------------------------------------------------------
  640 IF (KFLAG .EQ. -7) GO TO 670
      RH = .1
      RH = MAX(HMIN/ ABS(H),RH)
      H = H*RH
      IER = 0
      CALL DIFFUN (N, T, Y, SAVE1, IER, PW, IPIV, WORK, IWORK)
      IF ( IER .NE. 0 ) GO TO 685
      NJE = NJE + 1
      DO 650 I = 1,N
  650   Y(I,2) = H*SAVE1(I)
      IWEVAL = MITER
      IDOUB = 10
      IF (NQ .EQ. 1) GO TO 200
      NQ = 1
      L = 2
      IRET = 3
      GO TO 130
C-----------------------------------------------------------------------
C ALL RETURNS ARE MADE THROUGH THIS SECTION.  H IS SAVED IN HOLD
C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP.
C-----------------------------------------------------------------------
  660 KFLAG = -1
      GO TO 700
  670 KFLAG = -2
      GO TO 700
  680 KFLAG = -3
      GO TO 700
  685 KFLAG = -4
      GO TO 700
  690 RMAX = 10.
  700 HOLD = H
      JSTART = NQ
      RETURN
      END
        SUBROUTINE CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,
     *             NCLBLK,NBLOKS,BOTBLK,NRWBOT,PIVOT,IFLAG)
C
C***************************************************************
C
C  C R D C M P DECOMPOSES THE ALMOST BLOCK DIAGONAL MATRIX A
C  USING MODIFIED ALTERNATE ROW AND COLUMN ELIMINATION WITH
C  PARTIAL PIVOTING.  THE MATRIX  A  IS STORED IN THE ARRAYS
C  TOPBLK, ARRAY, AND BOTBLK.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C               N      - INTEGER
C                         THE ORDER OF THE LINEAR SYSTEM,
C                         GIVEN BY NBLOKS*NRWBLK + NOVRLP
C
C               TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C                         THE FIRST BLOCK OF THE ALMOST BLOCK
C                         DIAGONAL MATRIX A TO BE DECOMPOSED
C
C               NRWTOP - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C               NOVRLP - INTEGER
C                         THE NUMBER OF COLUMNS IN WHICH SUCC-
C                         ESSIVE BLOCKS OVERLAP, WHERE
C                                NOVRLP = NRWTOP + NRWBOT
C
C               ARRAY  - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C                         ARRAY(,,K) CONTAINS THE K-TH NRWBLK
C                         BY NCLBLK BLOCK OF THE MATRIX A
C
C               NRWBLK - INTEGER
C                         NUMBER OF ROWS IN K-TH BLOCK
C
C               NCLBLK - INTEGER
C                         NUMBER OF COLUMNS IN K-TH BLOCK
C
C               NBLOKS - INTEGER
C                         NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C                         THE MATRIX A
C
C               BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C                         THE LAST BLOCK OF THE MATRIX A
C
C               NRWBOT - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C                PIVOT - INTEGER(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C               TOPBLK,ARRAY,BOTBLK - ARRAYS CONTAINING THE
C                        DESIRED DECOMPOSITION OF THE MATRIX A
C                        (IF IFLAG = 0)
C
C                PIVOT - INTEGER(N)
C                         RECORDS THE PIVOTING INDICES DETER-
C                         MINED IN THE DECOMPOSITION
C
C               IFLAG  - INTEGER
C                         =  1, IF INPUT PARAMETERS ARE INVALID
C                         = -1, IF MATRIX IS SINGULAR
C                         =  0, OTHERWISE
C
C***************************************************************
C
        DOUBLE PRECISION TOPBLK,ARRAY,BOTBLK
        DOUBLE PRECISION ROWMAX,ROWPIV,ROWMLT,COLMAX,COLPIV
        DOUBLE PRECISION SWAP,COLMLT,PIVMAX,ZERO,TEMPIV
        INTEGER PIVOT(1)
        DIMENSION TOPBLK(NRWTOP,1),ARRAY(NRWBLK,NCLBLK,1),
     *          BOTBLK(NRWBOT,1)
        DATA ZERO/0.0D0/
C
C***************************************************************
C
C          ****  DEFINE THE CONSTANTS USED THROUGHOUT  ****
C
C***************************************************************
C
        IFLAG = 0
        PIVMAX = ZERO
        NRWTP1 = NRWTOP+1
        NROWEL = NRWBLK-NRWTOP
        NRWEL1 = NROWEL+1
        NVRLP0 = NOVRLP-1
C
C***************************************************************
C
C          ****  CHECK VALIDITY OF THE INPUT PARAMETERS....
C
C               IF PARAMETERS ARE INVALID THEN TERMINATE AT 10;
C                                         ELSE CONTINUE AT 100.
C
C***************************************************************
C
        IF(N.NE.NBLOKS*NRWBLK+NOVRLP)GO TO 10
        IF(NOVRLP.NE.NRWTOP+NRWBOT)GO TO 10
        IF(NCLBLK.NE.NOVRLP+NRWBLK)GO TO 10
        IF(NOVRLP.GT.NRWBLK)GO TO 10
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          PARAMETERS ARE ACCEPTABLE - CONTINUE AT 100.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        GO TO 100
10      CONTINUE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          PARAMETERS ARE INVALID.  SET IFLAG = 1, AND TERMINATE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IFLAG = 1
        RETURN
100     CONTINUE
C
C***************************************************************
C
C               ****  FIRST, IN TOPBLK....
C
C***************************************************************
C
C          ***  APPLY NRWTOP COLUMN ELIMINATIONS WITH COLUMN
C                 PIVOTING ....
C
C***************************************************************
C
        DO 190 I = 1,NRWTOP
           IPLUS1 = I+1
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE COLUMN PIVOT AND PIVOT INDEX
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           IPVT = I
           COLMAX = ABS(TOPBLK(I,I))
           DO 110 J = IPLUS1,NOVRLP
              TEMPIV = ABS(TOPBLK(I,J))
              IF(TEMPIV.LE.COLMAX)GO TO 110
                 IPVT = J
                 COLMAX = TEMPIV
110        CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           IF(PIVMAX+COLMAX.EQ.PIVMAX)GO TO 1000
           PIVMAX = MAX(COLMAX,PIVMAX)
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE COLUMNS
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           PIVOT(I) = IPVT
           IF(IPVT.EQ.I)GO TO 140
              DO 120 L = I,NRWTOP
                 SWAP = TOPBLK(L,IPVT)
                 TOPBLK(L,IPVT) = TOPBLK(L,I)
                 TOPBLK(L,I) = SWAP
120           CONTINUE
              DO 130 L = 1,NRWBLK
                 SWAP = ARRAY(L,IPVT,1)
                 ARRAY(L,IPVT,1) = ARRAY(L,I,1)
                 ARRAY(L,I,1) = SWAP
130           CONTINUE
140        CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS AND PERFORM COLUMN
C                       ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           COLPIV = TOPBLK(I,I)
           DO 180 J = IPLUS1,NOVRLP
              COLMLT = TOPBLK(I,J)/COLPIV
              TOPBLK(I,J) = COLMLT
              IF(IPLUS1.GT.NRWTOP)GO TO 160
                 DO 150 L = IPLUS1,NRWTOP
                    TOPBLK(L,J) = TOPBLK(L,J)-COLMLT*TOPBLK(L,I)
150              CONTINUE
160           CONTINUE
              DO 170 L = 1,NRWBLK
                 ARRAY(L,J,1) = ARRAY(L,J,1)-COLMLT*ARRAY(L,I,1)
170           CONTINUE
180        CONTINUE
190     CONTINUE
C
C***************************************************************
C
C          ****  IN EACH BLOCK ARRAY(,,K)....
C
C***************************************************************
C
        INCR = 0
        DO 395 K = 1,NBLOKS
           KPLUS1 = K+1
C
C          *****************************************************
C
C          ***  FIRST APPLY NRWBLK-NRWTOP ROW ELIMINATIONS WITH
C                       ROW PIVOTING....
C
C          *****************************************************
C
           DO 270 J = NRWTP1,NRWBLK
              JPLUS1 = J+1
              JMINN = J-NRWTOP
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE ROW PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              IPVT = JMINN
              ROWMAX = ABS(ARRAY(JMINN,J,K))
              LOOP = JMINN+1
              DO 210 I = LOOP,NRWBLK
                 TEMPIV = ABS(ARRAY(I,J,K))
                 IF(TEMPIV.LE.ROWMAX)GO TO 210
                 IPVT = I
                 ROWMAX = TEMPIV
210           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              IF(PIVMAX+ROWMAX.EQ.PIVMAX)GO TO  1000
              PIVMAX = MAX(ROWMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE ROWS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              INCRJ = INCR+J
              PIVOT(INCRJ) = INCR+IPVT+NRWTOP
              IF(IPVT.EQ.JMINN)GO TO 230
                 DO 220 L = J,NCLBLK
                    SWAP = ARRAY(IPVT,L,K)
                    ARRAY(IPVT,L,K) = ARRAY(JMINN,L,K)
                    ARRAY(JMINN,L,K) = SWAP
220              CONTINUE
230           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLERS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              ROWPIV = ARRAY(JMINN,J,K)
              DO 240 I = LOOP,NRWBLK
                 ARRAY(I,J,K) = ARRAY(I,J,K)/ROWPIV
240           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               PERFORM ROW ELIMINATION WITH COLUMN INDEXING
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              DO 260 L = JPLUS1,NCLBLK
                 ROWMLT = ARRAY(JMINN,L,K)
                 DO 250 I = LOOP,NRWBLK
                    ARRAY(I,L,K) = ARRAY(I,L,K)
     *                                -ROWMLT*ARRAY(I,J,K)
250              CONTINUE
260           CONTINUE
270        CONTINUE
C
C          *****************************************************
C
C          ***  NOW APPLY NRWTOP COLUMN ELIMINATIONS WITH
C                      COLUMN PIVOTING....
C
C          *****************************************************
C
           DO 390 I = NRWEL1,NRWBLK
              IPLUSN = I+NRWTOP
              IPLUS1 = I+1
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE COLUMN PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              IPVT = IPLUSN
              COLMAX = ABS(ARRAY(I,IPVT,K))
              LOOP = IPLUSN+1
              DO 310 J = LOOP,NCLBLK
                 TEMPIV = ABS(ARRAY(I,J,K))
                 IF(TEMPIV.LE.COLMAX)GO TO 310
                 IPVT = J
                 COLMAX = TEMPIV
310           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              IF(PIVMAX+COLMAX.EQ.PIVMAX)GO TO 1000
              PIVMAX = MAX(COLMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE COLUMNS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              INCRN = INCR+IPLUSN
              PIVOT(INCRN) = INCR+IPVT
              IRWBLK = IPLUSN-NRWBLK
              IF(IPVT.EQ.IPLUSN)GO TO 340
                 DO 315 L = I,NRWBLK
                    SWAP = ARRAY(L,IPVT,K)
                    ARRAY(L,IPVT,K) = ARRAY(L,IPLUSN,K)
                    ARRAY(L,IPLUSN,K) = SWAP
315              CONTINUE
                 IPVBLK = IPVT-NRWBLK
                 IF(K.EQ.NBLOKS)GO TO 330
                    DO 320 L = 1,NRWBLK
                       SWAP = ARRAY(L,IPVBLK,KPLUS1)
                       ARRAY(L,IPVBLK,KPLUS1)
     *                                 = ARRAY(L,IRWBLK,KPLUS1)
                       ARRAY(L,IRWBLK,KPLUS1) = SWAP
320                 CONTINUE
                    GO TO 340
330              CONTINUE
                 DO 335 L = 1,NRWBOT
                    SWAP = BOTBLK(L,IPVBLK)
                    BOTBLK(L,IPVBLK) = BOTBLK(L,IRWBLK)
                    BOTBLK(L,IRWBLK) = SWAP
335              CONTINUE
340           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS AND PERFORM COLUMN
C                       ELIMINATION
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              COLPIV = ARRAY(I,IPLUSN,K)
              DO 380 J = LOOP,NCLBLK
                 COLMLT = ARRAY(I,J,K)/COLPIV
                 ARRAY(I,J,K) = COLMLT
                 IF(I.EQ.NRWBLK)GO TO 350
                    DO 345 L = IPLUS1,NRWBLK
                       ARRAY(L,J,K) = ARRAY(L,J,K)
     *                                -COLMLT*ARRAY(L,IPLUSN,K)
345                 CONTINUE
350              CONTINUE
                 JRWBLK = J-NRWBLK
                 IF(K.EQ.NBLOKS)GO TO 370
                    DO 360 L = 1,NRWBLK
                       ARRAY(L,JRWBLK,KPLUS1) =
     *                                  ARRAY(L,JRWBLK,KPLUS1)
     *                         -COLMLT*ARRAY(L,IRWBLK,KPLUS1)
360                 CONTINUE
                    GO TO 380
370              CONTINUE
                 DO 375 L = 1,NRWBOT
                    BOTBLK(L,JRWBLK) = BOTBLK(L,JRWBLK)
     *                              -COLMLT*BOTBLK(L,IRWBLK)
375              CONTINUE
380           CONTINUE
390        CONTINUE
           INCR = INCR + NRWBLK
395     CONTINUE
C
C***************************************************************
C
C          ****  FINALLY, IN BOTBLK....
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          ***  APPLY NRWBOT ROW ELIMINATIONS WITH ROW
C                  PIVOTING....
C
C               IF BOT HAS JUST ONE ROW GO TO 500
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF(NRWBOT.EQ.1)GO TO 500
           DO 470 J = NRWTP1,NVRLP0
              JPLUS1 = J+1
              JMINN = J-NRWTOP
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE ROW PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              IPVT = JMINN
              ROWMAX = ABS(BOTBLK(JMINN,J))
              LOOP = JMINN+1
              DO 410 I = LOOP,NRWBOT
                 TEMPIV = ABS(BOTBLK(I,J))
                 IF(TEMPIV.LE.ROWMAX) GO TO 410
                 IPVT = I
                 ROWMAX = TEMPIV
410           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              IF(PIVMAX+ROWMAX.EQ.PIVMAX)GO TO 1000
              PIVMAX = MAX(ROWMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE ROWS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              INCRJ = INCR+J
              PIVOT(INCRJ) = INCR+IPVT+NRWTOP
              IF(IPVT.EQ.JMINN)GO TO 430
                 DO 420 L = J,NOVRLP
                    SWAP = BOTBLK(IPVT,L)
                    BOTBLK(IPVT,L) = BOTBLK(JMINN,L)
                    BOTBLK(JMINN,L) = SWAP
420              CONTINUE
430           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              ROWPIV = BOTBLK(JMINN,J)
              DO 440 I = LOOP,NRWBOT
                 BOTBLK(I,J) = BOTBLK(I,J)/ROWPIV
440           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               PERFORM ROW ELIMINATION WITH COLUMN INDEXING
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              DO 460 L = JPLUS1,NOVRLP
                 ROWMLT = BOTBLK(JMINN,L)
                 DO 450 I = LOOP,NRWBOT
                    BOTBLK(I,L) = BOTBLK(I,L)-ROWMLT*BOTBLK(I,J)
450              CONTINUE
460           CONTINUE
470        CONTINUE
500     CONTINUE
C
C***************************************************************
C
C          DONE PROVIDED THE LAST ELEMENT IS NOT ZERO
C
C***************************************************************
C
        IF(PIVMAX+ABS(BOTBLK(NRWBOT,NOVRLP)).NE.PIVMAX) RETURN
C
C***************************************************************
C
C       ****  MATRIX IS SINGULAR - SET IFLAG = - 1.
C                                  TERMINATE AT 1000.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
1000    CONTINUE
        IFLAG = -1
        RETURN
        END
        SUBROUTINE CRSLVE(TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,
     *             NCLBLK,NBLOKS,BOTBLK,NRWBOT,PIVOT,B,X)
C
C***************************************************************
C
C  C R S L V E  SOLVES THE LINEAR SYSTEM
C                       A*X = B
C  USING THE DECOMPOSITION ALREADY GENERATED IN  C R D C M P.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY  ...
C
C               TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C                         OUTPUT FROM  C R D C M P
C
C               NOVRLP - INTEGER
C                         THE NUMBER OF COLUMNS IN WHICH SUCC-
C                         ESSIVE BLOCKS OVERLAP, WHERE
C                                NOVRLP = NRWTOP + NRWBOT
C
C               NRWTOP - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C               ARRAY  - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C                         OUTPUT FROM  C R D C M P
C
C               NRWBLK - INTEGER
C                         NUMBER OF ROWS IN K-TH BLOCK
C
C               NCLBLK - INTEGER
C                         NUMBER OF COLUMNS IN K-TH BLOCK
C
C               NBLOKS - INTEGER
C                         NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C                         THE MATRIX A
C
C               BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C                         OUTPUT FROM  C R D C M P
C
C               NRWBOT - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C                PIVOT - INTEGER(N)
C                         THE PIVOT VECTOR FROM  C R D C M P
C
C                    B - DOUBLE PRECISION(N)
C                         THE RIGHT HAND SIDE VECTOR
C
C                    X - DOUBLE PRECISION(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C                    X - DOUBLE PRECISION(N)
C                         THE SOLUTION VECTOR
C
C***************************************************************
C
        DOUBLE PRECISION TOPBLK,ARRAY,BOTBLK,X,B
        DOUBLE PRECISION DOTPRD,XJ,XINCRJ,BINCRJ,SWAP
        INTEGER PIVOT(1)
        DIMENSION TOPBLK(NRWTOP,1),ARRAY(NRWBLK,NCLBLK,1),
     *          BOTBLK(NRWBOT,1),B(1),X(1)
C
C***************************************************************
C
C          ****  DEFINE THE CONSTANTS USED THROUGHOUT  ****
C
C***************************************************************
C
*
        NRWTP1 = NRWTOP+1
        NRWBK1 = NRWBLK+1
        NVRLP1 = NOVRLP+1
        NRWTP0 = NRWTOP-1
        NRWBT1 = NRWBOT+1
        NROWEL = NRWBLK-NRWTOP
        NRWEL1 = NROWEL+1
        NVRLP0 = NOVRLP-1
        NBLKS1 = NBLOKS+1
        NBKTOP = NRWBLK+NRWTOP
C
C***************************************************************
C
C               ****  FORWARD RECURSION  ****
C
C***************************************************************
C
C          ***  FIRST, IN TOPBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD SOLUTION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 130 J = 1,NRWTOP
           X(J) = B(J)/TOPBLK(J,J)
           IF(J.EQ.NRWTOP)GO TO 120
              XJ = -X(J)
              LOOP = J+1
              DO 110 I = LOOP,NRWTOP
                 B(I) = B(I)+TOPBLK(I,J)*XJ
110           CONTINUE
120        CONTINUE
130     CONTINUE
C
C       ********************************************************
C
C          ***  IN EACH BLOCK ARRAY(,,K)....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        INCR = 0
        DO 280 K = 1,NBLOKS
           INCRTP = INCR+NRWTOP
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD MODIFICATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           DO 220 J = 1,NRWTOP
              INCRJ = INCR+J
              XINCRJ = -X(INCRJ)
              DO 210 I = 1,NRWBLK
                 B(INCRTP+I) = B(INCRTP+I)+ARRAY(I,J,K)*XINCRJ
210           CONTINUE
220        CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           DO 240 J = NRWTP1,NRWBLK
              INCRJ = INCR+J
              JPIVOT = PIVOT(INCRJ)
              IF(JPIVOT.EQ.INCRJ)GO TO 225
                 SWAP = B(INCRJ)
                 B(INCRJ) = B(JPIVOT)
                 B(JPIVOT) = SWAP
225           CONTINUE
              BINCRJ = -B(INCRJ)
              LOOP = J-NRWTP0
              DO 230 I = LOOP,NRWBLK
                 B(INCRTP+I) = B(INCRTP+I)+ARRAY(I,J,K)*BINCRJ
230           CONTINUE
240        CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD SOLUTION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           DO 270 J = NRWBK1,NBKTOP
              INCRJ = INCR+J
              JRWTOP = J -NRWTOP
              X(INCRJ) = B(INCRJ)/ARRAY(JRWTOP,J,K)
              IF(J.EQ.NBKTOP)GO TO 260
                 XINCRJ = -X(INCRJ)
                 LOOP = J-NRWTP0
                 DO 250 I = LOOP,NRWBLK
                    B(INCRTP+I) = B(INCRTP+I)+ARRAY(I,J,K)*XINCRJ
250              CONTINUE
260           CONTINUE
270        CONTINUE
           INCR = INCR+NRWBLK
280     CONTINUE
C
C       ********************************************************
C
C          ***  FINALLY, IN BOTBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD MODIFICATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        INCRTP = INCR+NRWTOP
        DO 320 J = 1,NRWTOP
           INCRJ = INCR+J
           XINCRJ = -X(INCRJ)
           DO 310 I = 1,NRWBOT
              B(INCRTP+I) = B(INCRTP+I)+BOTBLK(I,J)*XINCRJ
310        CONTINUE
320     CONTINUE
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD ELIMINATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF(NRWBOT.EQ.1)GO TO 350
           DO 340 J = NRWTP1,NVRLP0
              INCRJ = INCR+J
              JPIVOT = PIVOT(INCRJ)
              IF(JPIVOT.EQ.INCRJ)GO TO 325
                 SWAP = B(INCRJ)
                 B(INCRJ) = B(JPIVOT)
                 B(JPIVOT) = SWAP
325           CONTINUE
              BINCRJ = -B(INCRJ)
              LOOP = J-NRWTP0
              DO 330 I = LOOP,NRWBOT
                 B(INCRTP+I) = B(INCRTP+I)+BOTBLK(I,J)*BINCRJ
330           CONTINUE
340        CONTINUE
350     CONTINUE
C
C***************************************************************
C
C               ****  BACKWARD RECURSION  ****
C
C***************************************************************
C
C          ***  FIRST IN BOTBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD SOLUTION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 430 LL = 1,NRWBOT
           J = NVRLP1-LL
           INCRJ = INCR+J
           NRWBTL = NRWBT1-LL
           X(INCRJ) = B(INCRJ)/BOTBLK(NRWBTL,J)
           IF(LL.EQ.NRWBOT)GO TO 420
              XINCRJ = -X(INCRJ)
              LOOP = NRWBOT-LL
              DO 410 I = 1,LOOP
                 B(INCRTP+I) = B(INCRTP+I)+BOTBLK(I,J)*XINCRJ
410           CONTINUE
420        CONTINUE
430     CONTINUE
C
C       ********************************************************
C
C          ***  THEN IN EACH BLOCK ARRAY(,,K)....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 490 L = 1,NBLOKS
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           K = NBLKS1-L
           INCR = INCR-NRWBLK
           DO 450 L1 = NRWEL1,NRWBLK
              I = NRWBLK+NRWEL1-L1
              IPLUSN = I+NRWTOP
              LOOP = IPLUSN+1
              INCRN = INCR+IPLUSN
              DOTPRD = X(INCRN)
              DO 440 J = LOOP,NCLBLK
                 INCRJ = INCR+J
                 DOTPRD = DOTPRD-ARRAY(I,J,K)*X(INCRJ)
440           CONTINUE
              X(INCRN) = DOTPRD
              IPVTN = PIVOT(INCRN)
              IF(INCRN.EQ.IPVTN)GO TO 445
                 SWAP = X(INCRN)
                 X(INCRN) = X(IPVTN)
                 X(IPVTN) = SWAP
445           CONTINUE
450        CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD MODIFICATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           INCRTP = INCR+NRWTOP
           DO 460 J = NRWBK1,NCLBLK
              INCRJ = INCR+J
              XINCRJ = -X(INCRJ)
              DO 455 I = 1,NROWEL
                 B(INCRTP+I) = B(INCRTP+I)+ARRAY(I,J,K)*XINCRJ
455           CONTINUE
460        CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD SOLUTION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           DO 480 LL = 1,NROWEL
              J = NRWBK1-LL
              INCRJ = INCR+J
              NRWELL = NRWEL1-LL
              X(INCRJ) = B(INCRJ)/ARRAY(NRWELL,J,K)
              IF(LL.EQ.NROWEL)GO TO 470
                 XINCRJ = -X(INCRJ)
                 LOOP = NROWEL-LL
                 DO 465 I = 1,LOOP
                    B(INCRTP+I) = B(INCRTP+I)+ARRAY(I,J,K)*XINCRJ
465              CONTINUE
470           CONTINUE
480        CONTINUE
490     CONTINUE
C
C       ********************************************************
C
C          ***  IN TOPBLK FINISH WITH....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD ELIMINATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 520 L = 1,NRWTOP
           I = NRWTP1-L
           LOOP = I+1
           DOTPRD = X(I)
           DO 510 J = LOOP,NOVRLP
              DOTPRD = DOTPRD-TOPBLK(I,J)*X(J)
510        CONTINUE
           X(I) = DOTPRD
           IPVTI = PIVOT(I)
           IF(I.EQ.IPVTI)GO TO 515
                 SWAP = X(I)
                 X(I) = X(IPVTI)
                 X(IPVTI) = SWAP
515        CONTINUE
520     CONTINUE
        RETURN
        END
      BLOCK DATA
      IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C IN THE FOLLOWING DATA STATEMENT, SET..
C    LOUT   = THE LOGICAL UNIT NUMBER FOR THE OUTPUT OF MESSAGES DURING
C             THE INTEGRATION.
C    NOGAUS = SET .EQ. 1 IF THE GAUSS-LEGENDRE COLLOCATION POINTS ARE
C             NOT DESIRED WHEN NCC = 2 (SEE PDECOL AND COLPNT).
C    MAXDER = SET .EQ. 5.  ITS VALUE REPRESENTS THE MAXIMUM ORDER OF
C             THE TIME INTEGRATION ALLOWED.  ITS VALUE AFFECTS THE STOR-
C             AGE REQUIRED IN /HISTRY/ AND MAY BE CHANGED IF DESIRED
C             (SEE COSET FOR RESTRICTIONS).
C    UROUND = THE UNIT ROUNDOFF OF THE MACHINE, I.E. THE SMALLEST
C             POSITIVE U SUCH THAT 1. + U .NE. 1. ON THE MACHINE.
C-----------------------------------------------------------------------
      COMMON /GEAR1/ DUM(5),UROUND,IDUM(4)
      COMMON /OPTION/ NOGAUS,MAXDER
      COMMON /IOUNIT/ LOUT
      DATA LOUT,NOGAUS,MAXDER,UROUND/6,0,5,2.22D-16/
      END
      SUBROUTINE BSPLVD ( XT, K, X, ILEFT, VNIKX, NDERIV )
      IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C THIS SUBROUTINE IS PART OF THE B-SPLINE PACKAGE FOR THE STABLE
C EVALUATION OF ANY B-SPLINE BASIS FUNCTION OR DERIVATIVE VALUE.
C SEE REFERENCE BELOW.
C
C CALCULATES THE VALUE AND THE FIRST NDERIV-1 DERIVATIVES OF ALL
C B-SPLINES WHICH DO NOT VANISH AT X.  THE ROUTINE FILLS THE TWO-
C DIMENSIONAL ARRAY VNIKX(J,IDERIV), J=IDERIV, ... ,K WITH NONZERO
C VALUES OF B-SPLINES OF ORDER K+1-IDERIV, IDERIV=NDERIV, ... ,1, BY
C REPEATED CALLS TO BSPLVN.
C
C REFERENCE
C
C    DEBOOR, C., PACKAGE FOR CALCULATING WITH B-SPLINES, SIAM J.
C      NUMER. ANAL., VOL. 14, NO. 3, JUNE 1977, PP. 441-472.
C
C PACKAGE ROUTINES CALLED..  BSPLVN
C USER ROUTINES CALLED..     NONE
C CALLED BY..                COLPNT,INITAL,VALUES
C FORTRAN FUNCTIONS USED..   FLOAT,MAX0
C-----------------------------------------------------------------------
      DIMENSION XT(1),VNIKX(K,NDERIV)
      DIMENSION A(20,20)
      KO = K + 1 - NDERIV
      CALL BSPLVN(XT,KO,1,X,ILEFT,VNIKX(NDERIV,NDERIV))
      IF (NDERIV .LE. 1) GO TO 120
      IDERIV = NDERIV
      DO 20 I=2,NDERIV
        IDERVM = IDERIV-1
        DO 10 J=IDERIV,K
   10     VNIKX(J-1,IDERVM) = VNIKX(J,IDERIV)
        IDERIV = IDERVM
        CALL BSPLVN(XT,0,2,X,ILEFT,VNIKX(IDERIV,IDERIV))
   20   CONTINUE
      DO 40 I=1,K
        DO 30 J=1,K
   30     A(I,J) = 0.
   40   A(I,I) = 1.
      KMD = K
      DO 110 M=2,NDERIV
        KMD = KMD - 1
        FKMD =  FLOAT(KMD)
        I = ILEFT
        J = K
   50   JM1 = J-1
        IPKMD = I + KMD
        DIFF = XT(IPKMD) -XT(I)
        IF (JM1 .EQ. 0) GO TO 80
        IF (DIFF .EQ. 0.) GO TO 70
        DO 60 L=1,J
   60     A(L,J) = (A(L,J) - A(L,J-1))/DIFF*FKMD
   70   J = JM1
        I = I - 1
        GO TO 50
   80   IF (DIFF .EQ. 0.) GO TO 90
        A(1,1) = A(1,1)/DIFF*FKMD
   90   DO 110 I=1,K
          V = 0.
          JLOW = MAX0(I,M)
          DO 100 J=JLOW,K
  100       V = A(I,J)*VNIKX(J,M) + V
  110     VNIKX(I,M) = V
  120 RETURN
      END
      SUBROUTINE BSPLVN ( XT, JHIGH, INDEX, X, ILEFT, VNIKX )
      IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C THIS SUBROUTINE IS PART OF THE B-SPLINE PACKAGE FOR THE STABLE
C EVALUATION OF ANY B-SPLINE BASIS FUNCTION OR DERIVATIVE VALUE.
C SEE REFERENCE BELOW.
C
C CALCULATES THE VALUE OF ALL POSSIBLY NONZERO B-SPLINES AT THE
C POINT X OF ORDER MAX(JHIGH,(J+1)(INDEX-1)) FOR THE BREAKPOINT SEQ-
C UENCE XT.  ASSUMING THAT XT(ILEFT) .LE. X .LE. XT(ILEFT+1), THE ROUT-
C INE RETURNS THE B-SPLINE VALUES IN THE ONE DIMENSIONAL ARRAY VNIKX.
C
C REFERENCE
C
C    DEBOOR, C., PACKAGE FOR CALCULATING WITH B-SPLINES, SIAM J.
C      NUMER. ANAL., VOL. 14, NO. 3, JUNE 1977, PP. 441-472.
C
C PACKAGE ROUTINES CALLED..  NONE
C USER ROUTINES CALLED..     NONE
C CALLED BY..                BSPLVD
C FORTRAN FUNCTIONS USED..   NONE
C-----------------------------------------------------------------------
      DIMENSION XT(1),VNIKX(1)
      DIMENSION DELTAM(20),DELTAP(20)
      DATA J/1/,DELTAM/20*0.D+0/,DELTAP/20*0.D+0/
      GO TO (10,20),INDEX
   10 J = 1
      VNIKX(1) = 1.
      IF (J .GE. JHIGH) GO TO 40
   20 IPJ = ILEFT+J
      DELTAP(J) = XT(IPJ) - X
      IMJP1 = ILEFT-J+1
      DELTAM(J) = X - XT(IMJP1)
      VMPREV = 0.
      JP1 = J+1
      DO 30 L=1,J
        JP1ML = JP1-L
        VM = VNIKX(L)/(DELTAP(L) + DELTAM(JP1ML))
        VNIKX(L) = VM*DELTAP(L) + VMPREV
   30   VMPREV = VM*DELTAM(JP1ML)
      VNIKX(JP1) = VMPREV
      J = JP1
      IF (J .LT. JHIGH) GO TO 20
   40 RETURN
      END
      SUBROUTINE INTERV ( XT, LXT, X, ILEFT, MFLAG )
      IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C THIS SUBROUTINE IS PART OF THE B-SPLINE PACKAGE FOR THE STABLE
C EVALUATION OF ANY B-SPLINE BASIS FUNCTION OR DERIVATIVE VALUE.
C SEE REFERENCE BELOW.
C
C COMPUTES LARGEST ILEFT IN (1,LXT) SUCH THAT XT(ILEFT) .LE. X.  THE
C PROGRAM STARTS THE SEARCH FOR ILEFT WITH THE VALUE OF ILEFT THAT WAS
C RETURNED AT THE PREVIOUS CALL (AND WAS SAVED IN THE LOCAL VARIABLE
C ILO) TO MINIMIZE THE WORK IN THE COMMON CASE THAT THE VALUE OF X ON
C THIS CALL IS CLOSE TO THE VALUE OF X ON THE PREVIOUS CALL.  SHOULD
C THIS ASSUMPTION NOT BE VALID, THEN THE PROGRAM LOCATES ILO AND IHI
C SUCH THAT XT(ILO) .LE. X .LT. XT(IHI) AND, ONCE THEY ARE FOUND USES
C BISECTION TO FIND THE CORRECT VALUE FOR ILEFT.
C
C REFERENCE
C
C    DEBOOR, C., PACKAGE FOR CALCULATING WITH B-SPLINES, SIAM J.
C      NUMER. ANAL., VOL. 14, NO. 3, JUNE 1977, PP. 441-472.
C
C PACKAGE ROUTINES CALLED..  NONE
C USER ROUTINES CALLED..     NONE
C CALLED BY..                COLPNT,INITAL,VALUES
C FORTRAN FUNCTIONS USED..   NONE
C-----------------------------------------------------------------------
      DIMENSION XT(LXT)
      IF(MFLAG.EQ.-2) ILO = 1
      IHI = ILO + 1
      IF (IHI .LT. LXT) GO TO 20
      IF (X .GE. XT(LXT)) GO TO 110
      IF (LXT .LE. 1) GO TO 90
      ILO = LXT - 1
      GO TO 21
   20 IF (X .GE. XT(IHI)) GO TO 40
   21 IF (X .GE. XT(ILO)) GO TO 100
C-----------------------------------------------------------------------
C NOW X .LT. XT(IHI).  FIND LOWER BOUND.
C-----------------------------------------------------------------------
   30 ISTEP = 1
   31 IHI = ILO
      ILO = IHI - ISTEP
      IF (ILO .LE. 1) GO TO 35
      IF (X .GE. XT(ILO)) GO TO 50
      ISTEP = ISTEP*2
      GO TO 31
   35 ILO = 1
      IF (X .LT. XT(1)) GO TO 90
      GO TO 50
C-----------------------------------------------------------------------
C NOW X .GE. XT(ILO).  FIND UPPER BOUND.
C-----------------------------------------------------------------------
   40 ISTEP = 1
   41 ILO = IHI
      IHI = ILO + ISTEP
      IF (IHI .GE. LXT) GO TO 45
      IF (X .LT. XT(IHI)) GO TO 50
      ISTEP = ISTEP*2
      GO TO 41
   45 IF (X .GE. XT(LXT)) GO TO 110
      IHI = LXT
C-----------------------------------------------------------------------
C NOW XT(ILO) .LE. X .LT. XT(IHI).  NARROW THE INTERVAL.
C-----------------------------------------------------------------------
   50 MIDDLE = (ILO + IHI)/2
      IF (MIDDLE .EQ. ILO) GO TO 100
C-----------------------------------------------------------------------
C NOTE..  IT IS ASSUMED THAT MIDDLE = ILO IN CASE IHI = ILO+1.
C-----------------------------------------------------------------------
      IF (X .LT. XT(MIDDLE)) GO TO 53
      ILO = MIDDLE
      GO TO 50
   53 IHI = MIDDLE
      GO TO 50
C-----------------------------------------------------------------------
C SET OUTPUT AND RETURN.
C-----------------------------------------------------------------------
   90 MFLAG = -1
      ILEFT = 1
      RETURN
  100 MFLAG = 0
      ILEFT = ILO
      RETURN
  110 MFLAG = 1
      ILEFT = LXT
      RETURN
      END
      SUBROUTINE COLPNT(X, XC, XT)
      IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C COLPNT IS CALLED ONLY ONCE BY INITAL TO DEFINE THE REQUIRED COLLOCA-
C TION POINTS WHICH ARE TO BE USED WITH THE USER SELECTED PIECEWISE
C POLYNOMIAL SPACE.  THE COLLOCATION POINTS ARE CHOSEN SUCH THAT THEY
C ARE EITHER THE POINTS AT WHICH THE PIECEWISE POLYNOMIAL SPACE BASIS
C FUNCTIONS ATTAIN THEIR UNIQUE MAXIMUM VALUES, OR, THE GAUSS-LEGENDRE
C QUADRATURE POINTS WITHIN EACH PIECEWISE POLYNOMIAL SPACE SUBINTERVAL,
C DEPENDING UPON THE SPACE BEING USED AND THE DESIRE OF THE USER.
C
C PACKAGE ROUTINES CALLED..  BSPLVD,INTERV
C USER ROUTINES CALLED..     NONE
C CALLED BY..                INITAL
C FORTRAN FUNCTIONS USED..   NONE
C-----------------------------------------------------------------------
      DIMENSION RHO(40),X(1),XC(1),XT(1)
      COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD
      COMMON /OPTION/ NOGAUS,MAXDER
      DATA ILEFT/0/
C-----------------------------------------------------------------------
C IF THE VARIABLE NOGAUS IN THE COMMON BLOCK /OPTION/ IS SET .EQ. 1,
C THE USE OF THE GAUSS-LEGENDRE POINTS IS PROHIBITED FOR ALL CASES.
C NOGAUS IS CURRENTLY SET .EQ. 0 BY A DATA STATEMENT IN THE BLOCK DATA.
C THE USER MAY CHANGE THIS AS DESIRED.
C-----------------------------------------------------------------------
      IF ( NCC .NE. 2  .OR. NOGAUS .EQ. 1 ) GO TO 200
C-----------------------------------------------------------------------
C COMPUTE THE COLLOCATION POINTS TO BE AT THE GAUSS-LEGENDRE POINTS IN
C EACH PIECEWISE POLYNOMIAL SPACE SUBINTERVAL.  THE ARRAY RHO IS SET TO
C CONTAIN THE GAUSS-LEGENDRE POINTS FOR THE STANDARD INTERVAL (-1,1).
C-----------------------------------------------------------------------
      IPTS = KORD - 2
      GO TO (10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,
     *    180),IPTS
   10 RHO(1) = 0.
      GO TO 190
   20 RHO(2) = .577350269189626D-00
      RHO(1) = - RHO(2)
      GO TO 190
   30 RHO(3) = .774596669241483D-00
      RHO(1) = - RHO(3)
      RHO(2) = 0.
      GO TO 190
   40 RHO(3) = .339981043584856D-00
      RHO(2) = - RHO(3)
      RHO(4) = .861136311594053D-00
      RHO(1) = - RHO(4)
      GO TO 190
   50 RHO(4) = .538469310105683D-00
      RHO(2) = - RHO(4)
      RHO(5) = .906179845938664D-00
      RHO(1) = - RHO(5)
      RHO(3) = 0.
      GO TO 190
   60 RHO(4) = .238619186083197D-00
      RHO(3) = - RHO(4)
      RHO(5) = .661209386466265D-00
      RHO(2) = - RHO(5)
      RHO(6) = .932469514203152D-00
      RHO(1) = - RHO(6)
      GO TO 190
   70 RHO(5) = .405845151377397D-00
      RHO(3) = - RHO(5)
      RHO(6) = .741531185599394D-00
      RHO(2) = - RHO(6)
      RHO(7) = .949107912342759D-00
      RHO(1) = - RHO(7)
      RHO(4) = 0.
      GO TO 190
   80 RHO(5) = .183434642495650D-00
      RHO(4) = - RHO(5)
      RHO(6) = .525532409916329D-00
      RHO(3) = - RHO(6)
      RHO(7) = .796666477413627D-00
      RHO(2) = - RHO(7)
      RHO(8) = .960289856497536D-00
      RHO(1) = - RHO(8)
      GO TO 190
   90 RHO( 5) = .0
      RHO( 6) = .324253423403809D-00
      RHO( 7) = .613371432700590D-00
      RHO( 8) = .836031107326636D-00
      RHO( 9) = .968160239507626D-00
      DO 95 I=1,4
   95   RHO(I) = -RHO(10-I)
      GO TO 190
  100 RHO( 6) = .148874338981631D-00
      RHO( 7) = .433395394129247D-00
      RHO( 8) = .679409568299024D-00
      RHO( 9) = .865063366688984D-00
      RHO(10) = .973906528517172D-00
      DO 105 I=1,5
  105   RHO(I) = -RHO(11-I)
      GO TO 190
  110 RHO( 6) = .0
      RHO( 7) = .269543155952345D-00
      RHO( 8) = .519096129206812D-00
      RHO( 9) = .730152005574049D-00
      RHO(10) = .887062599768095D-00
      RHO(11) = .978228658146057D-00
      DO 115 I=1,5
  115   RHO(I) = -RHO(12-I)
      GO TO 190
  120 RHO( 7) = .125233408511469D-00
      RHO( 8) = .367831498998180D-00
      RHO( 9) = .587317954286617D-00
      RHO(10) = .769902674194305D-00
      RHO(11) = .904117256370475D-00
      RHO(12) = .981560634246719D-00
      DO 125 I=1,6
  125   RHO(I) = -RHO(13-I)
      GO TO 190
  130 RHO( 7) = .0
      RHO( 8) = .230458315955135D-00
      RHO( 9) = .448492751036447D-00
      RHO(10) = .642349339440340D-00
      RHO(11) = .801578090733310D-00
      RHO(12) = .917598399222978D-00
      RHO(13) = .984183054718588D-00
      DO 135 I=1,6
  135   RHO(I) = -RHO(14-I)
      GO TO 190
  140 RHO( 8) = .108054948707344D-00
      RHO( 9) = .319112368927890D-00
      RHO(10) = .515248636358154D-00
      RHO(11) = .687292904811685D-00
      RHO(12) = .827201315069765D-00
      RHO(13) = .928434883663574D-00
      RHO(14) = .986283808696812D-00
      DO 145 I=1,7
  145   RHO(I) = -RHO(15-I)
      GO TO 190
  150 RHO( 8) = .0
      RHO( 9) = .201194093997435D-00
      RHO(10) = .394151347077563D-00
      RHO(11) = .570972172608539D-00
      RHO(12) = .724417731360170D-00
      RHO(13) = .848206583410427D-00
      RHO(14) = .937273392400706D-00
      RHO(15) = .987992518020485D-00
      DO 155 I = 1,7
  155   RHO(I) = -RHO(16-I)
      GO TO 190
  160 RHO( 9) = .950125098376374D-01
      RHO(10) = .281603550779259D-00
      RHO(11) = .458016777657227D-00
      RHO(12) = .617876244402644D-00
      RHO(13) = .755404408355003D-00
      RHO(14) = .865631202387832D-00
      RHO(15) = .944575023073233D-00
      RHO(16) = .989400934991650D-00
      DO 165 I=1,8
  165   RHO(I) = -RHO(17-I)
      GO TO 190
  170 RHO( 9) = .0
      RHO(10) = .178484181495848D-00
      RHO(11) = .351231763453876D-00
      RHO(12) = .512690537086477D-00
      RHO(13) = .657671159216691D-00
      RHO(14) = .781514003896801D-00
      RHO(15) = .880239153726986D-00
      RHO(16) = .950675521768768D-00
      RHO(17) = .990575475314417D-00
      DO 175 I=1,8
  175   RHO(I) = -RHO(18-I)
      GO TO 190
  180 RHO(10) = .847750130417353D-01
      RHO(11) = .251886225691506D-00
      RHO(12) = .411751161462843D-00
      RHO(13) = .559770831073948D-00
      RHO(14) = .691687043060353D-00
      RHO(15) = .803704958972523D-00
      RHO(16) = .892602466497556D-00
      RHO(17) = .955823949571398D-00
      RHO(18) = .991565168420931D-00
      DO 185 I=1,9
  185   RHO(I) = -RHO(19-I)
C-----------------------------------------------------------------------
C COMPUTE THE GAUSS-LEGENDRE COLLOCATION POINTS IN EACH SUBINTERVAL.
C-----------------------------------------------------------------------
  190 DO 195 I=1,NINT
        FAC = ( X(I+1) - X(I) ) * .5
        DO 195 J = 1,IPTS
          KNOT = IPTS * (I-1) + J + 1
  195     XC(KNOT) = X(I) + FAC * ( RHO(J) + 1. )
      XC(1) = X(1)
      XC(NCPTS) = X(NINT+1)
      RETURN
C-----------------------------------------------------------------------
C COMPUTE THE COLLOCATION POINTS TO BE AT THE POINTS WHERE THE BASIS
C FUNCTIONS ATTAIN THEIR MAXIMA.  A BISECTION METHOD IS USED TO FIND
C THE POINTS TO MACHINE PRECISION.  THIS PROCESS COULD BE SPEEDED UP
C BY USING A SECANT METHOD IF DESIRED.
C-----------------------------------------------------------------------
  200 ITOP = NCPTS - 1
      MFLAG = -2
      XC(1) = X(1)
      XC(NCPTS) = X(NINT+1)
      DO 240 I=2,ITOP
        XOLD = 1.D+20
        XL = XT(I)
        XR = XT(I+KORD)
  210   XNEW = .5 * (XL + XR)
        IF( XOLD .EQ. XNEW ) GO TO 240
        CALL INTERV(XT,NCPTS,XNEW,ILEFT,MFLAG)
        CALL BSPLVD(XT,KORD,XNEW,ILEFT,RHO,2)
        DO 220 J=1,KORD
          IF( I .EQ. J + ILEFT - KORD ) GO TO 230
  220   CONTINUE
  230   XVAL = RHO(KORD+J)
        IF( XVAL .EQ. 0.0 ) XR = XNEW
        IF( XVAL .GT. 0.0 ) XL = XNEW
        IF( XVAL .LT. 0.0 ) XR = XNEW
        XOLD = XNEW
        GO TO 210
  240 XC(I) = XR
      RETURN
      END
      SUBROUTINE COSET (METH, NQ, EL, TQ)
      IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C COSET IS CALLED BY THE INTEGRATOR AND SETS COEFFICIENTS USED THERE.
C THE VECTOR EL, OF LENGTH NQ + 1, DETERMINES THE BASIC METHOD.
C THE VECTOR TQ, OF LENGTH 4, IS INVOLVED IN ADJUSTING THE STEP SIZE
C IN RELATION TO TRUNCATION ERROR.  ITS VALUES ARE GIVEN BY THE
C PERTST ARRAY.
C
C THE VECTORS EL AND TQ DEPEND ON METH AND NQ.
C THE MAXIMUM ORDER, MAXDER, OF THE METHODS AVAILABLE IS CURRENTLY
C 12 FOR THE ADAMS METHODS AND 5 FOR THE BDF METHODS.  MAXDER DEFAULTS
C TO 5 UNLESS THE USER SETS MAXDER TO SOME OTHER LEGITIMATE VALUE
C THROUGH THE COMMON BLOCK /OPTION/.  SEE PDECOL FOR ADDITIONAL DETAILS.
C LMAX = MAXDER + 1 IS THE NUMBER OF COLUMNS IN THE Y ARRAY (SEE STIFIB
C AND THE VARIABLE C, Y, OR WORK(IW10) IN PDECOL).
C
C THE COEFFICIENTS IN PERTST NEED BE GIVEN TO ONLY ABOUT
C ONE PERCENT ACCURACY.  THE ORDER IN WHICH THE GROUPS APPEAR BELOW
C IS..  COEFFICIENTS FOR ORDER NQ - 1, COEFFICIENTS FOR ORDER NQ,
C COEFFICIENTS FOR ORDER NQ + 1.  WITHIN EACH GROUP ARE THE
C COEFFICIENTS FOR THE ADAMS METHODS, FOLLOWED BY THOSE FOR THE
C BDF METHODS.
C
C REFERENCE
C
C   GEAR, C.W., NUMERICAL INITIAL VALUE PROBLEMS IN ORDINARY
C     DIFFERENTIAL EQUATIONS, PRENTICE-HALL, ENGLEWOOD CLIFFS,
C     N. J., 1971.
C
C PACKAGE ROUTINES CALLED..  NONE
C USER ROUTINES CALLED..     NONE
C CALLED BY..                STIFIB
C FORTRAN FUNCTIONS USED..   FLOAT
C-----------------------------------------------------------------------
      DIMENSION PERTST(12,2,3),EL(13),TQ(4)
      DATA  PERTST / 1.,1.,2.,1.,.3158,.07407,.01391,.002182,
     *                 .0002945,.00003492,.000003692,.0000003524,
     *               1.,1.,.5,.1667,.04167,1.,1.,1.,1.,1.,1.,1.,
     *               2.,12.,24.,37.89,53.33,70.08,87.97,106.9,
     *                 126.7,147.4,168.8,191.0,
     *               2.0,4.5,7.333,10.42,13.7,1.,1.,1.,1.,1.,1.,1.,
     *               12.0,24.0,37.89,53.33,70.08,87.97,106.9,
     *                 126.7,147.4,168.8,191.0,1.,
     *               3.0,6.0,9.167,12.5,1.,1.,1.,1.,1.,1.,1.,1. /
C
      GO TO (1,2),METH
    1 GO TO (101,102,103,104,105,106,107,108,109,110,111,112),NQ
    2 GO TO (201,202,203,204,205),NQ
C-----------------------------------------------------------------------
C THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO MACHINE ACCURACY.
C FOR A GIVEN ORDER NQ, THEY CAN BE CALCULATED BY USE OF THE
C GENERATING POLYNOMIAL L(T), WHOSE COEFFICIENTS ARE EL(I)..
C      L(T) = EL(1) + EL(2)*T + ... + EL(NQ+1)*T**NQ.
C FOR THE IMPLICIT ADAMS METHODS, L(T) IS GIVEN BY
C      DL/DT = (T+1)*(T+2)* ... *(T+NQ-1)/K,    L(-1) = 0,
C WHERE                 K = FACTORIAL(NQ-1).
C FOR THE BDF METHODS,
C      L(T) = (T+1)*(T+2)* ... *(T+NQ)/K,
C WHERE         K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ).
C
C THE ORDER IN WHICH THE GROUPS APPEAR BELOW IS..
C IMPLICIT ADAMS METHODS OF ORDERS 1 TO 12,
C BDF METHODS OF ORDERS 1 TO 5.
C-----------------------------------------------------------------------
 101  EL(1) = 1.0D-00
      GO TO 900
 102  EL(1) = 0.5D-00
      EL(3) = 0.5D-00
      GO TO 900
 103  EL(1) = 4.1666666666667D-01
      EL(3) = 0.75D-00
      EL(4) = 1.6666666666667D-01
      GO TO 900
 104  EL(1) = 0.375D-00
      EL(3) = 9.1666666666667D-01
      EL(4) = 3.3333333333333D-01
      EL(5) = 4.1666666666667D-02
      GO TO 900
 105  EL(1) = 3.4861111111111D-01
      EL(3) = 1.0416666666667D-00
      EL(4) = 4.8611111111111D-01
      EL(5) = 1.0416666666667D-01
      EL(6) = 8.3333333333333D-03
      GO TO 900
 106  EL(1) = 3.2986111111111D-01
      EL(3) = 1.1416666666667D-00
      EL(4) = 0.625D-00
      EL(5) = 1.7708333333333D-01
      EL(6) = 0.025D-00
      EL(7) = 1.3888888888889D-03
      GO TO 900
 107  EL(1) = 3.1559193121693D-01
      EL(3) = 1.225D-00
      EL(4) = 7.5185185185185D-01
      EL(5) = 2.5520833333333D-01
      EL(6) = 4.8611111111111D-02
      EL(7) = 4.8611111111111D-03
      EL(8) = 1.9841269841270D-04
      GO TO 900
 108  EL(1) = 3.0422453703704D-01
      EL(3) = 1.2964285714286D-00
      EL(4) = 8.6851851851852D-01
      EL(5) = 3.3576388888889D-01
      EL(6) = 7.7777777777778D-02
      EL(7) = 1.0648148148148D-02
      EL(8) = 7.9365079365079D-04
      EL(9) = 2.4801587301587D-05
      GO TO 900
 109  EL(1) = 2.9486800044092D-01
      EL(3) = 1.3589285714286D-00
      EL(4) = 9.7655423280423D-01
      EL(5) = 0.4171875D-00
      EL(6) = 1.1135416666667D-01
      EL(7) = 0.01875D-00
      EL(8) = 1.9345238095238D-03
      EL(9) = 1.1160714285714D-04
      EL(10)= 2.7557319223986D-06
      GO TO 900
 110  EL(1) = 2.8697544642857D-01
      EL(3) = 1.4144841269841D-00
      EL(4) = 1.0772156084656D-00
      EL(5) = 4.9856701940035D-01
      EL(6) = 0.1484375D-00
      EL(7) = 2.9060570987654D-02
      EL(8) = 3.7202380952381D-03
      EL(9) = 2.9968584656085D-04
      EL(10)= 1.3778659611993D-05
      EL(11)= 2.7557319223986D-07
      GO TO 900
 111  EL(1) = 2.8018959644394D-01
      EL(3) = 1.4644841269841D-00
      EL(4) = 1.1715145502646D-00
      EL(5) = 5.7935819003527D-01
      EL(6) = 1.8832286155203D-01
      EL(7) = 4.1430362654321D-02
      EL(8) = 6.2111441798942D-03
      EL(9) = 6.2520667989418D-04
      EL(10)= 4.0417401528513D-05
      EL(11)= 1.5156525573192D-06
      EL(12)= 2.5052108385442D-08
      GO TO 900
 112  EL(1) = 2.7426554003160D-01
      EL(3) = 1.5099386724387D-00
      EL(4) = 1.2602711640212D-00
      EL(5) = 6.5923418209877D-01
      EL(6) = 2.3045800264550D-01
      EL(7) = 5.5697246105232D-02
      EL(8) = 9.4394841269841D-03
      EL(9) = 1.1192749669312D-03
      EL(10)= 9.0939153439153D-05
      EL(11)= 4.8225308641975D-06
      EL(12)= 1.5031265031265D-07
      EL(13)= 2.0876756987868D-09
      GO TO 900
 201  EL(1) = 1.0D-00
      GO TO 900
 202  EL(1) = 6.6666666666667D-01
      EL(3) = 3.3333333333333D-01
      GO TO 900
 203  EL(1) = 5.4545454545455D-01
      EL(3) = EL(1)
      EL(4) = 9.0909090909091D-02
      GO TO 900
 204  EL(1) = 0.48D-00
      EL(3) = 0.7D-00
      EL(4) = 0.2D-00
      EL(5) = 0.02D-00
      GO TO 900
 205  EL(1) = 4.3795620437956D-01
      EL(3) = 8.2116788321168D-01
      EL(4) = 3.1021897810219D-01
      EL(5) = 5.4744525547445D-02
      EL(6) = 3.6496350364964D-03
C
 900  DO 910 K = 1,3
 910    TQ(K) = PERTST(NQ,METH,K)
      TQ(4) = .5D-00*TQ(2)/ FLOAT(NQ+2)
      RETURN
      END
      SUBROUTINE DIFFF(T,X,IPT,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE,CMAX,
     *                 SAVE2)
      IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C SUBROUTINE DIFFF IS USED IF MITER=2 TO PROVIDE FINITE DIFFERENCE
C APPROXIMATIONS FOR THE PARTIAL DERIVATIVES OF THE K-TH USER DEFINED
C FUNCTION IN THE F ROUTINE WITH RESPECT TO THE VARIABLES U, UX, AND
C UXX.  THESE PARTIALS WITH RESPECT TO U, UX, AND UXX ARE COMPUTED,
C STORED, AND RETURNED IN THE NPDE BY NPDE ARRAYS DFDU, DFDUX, AND
C DFDUXX, RESPECTIVELY.
C
C PACKAGE ROUTINES CALLED..  NONE
C USER ROUTINES CALLED..     F
C CALLED BY..                PSETIB
C FORTRAN FUNCTIONS USED..   MAX
C-----------------------------------------------------------------------
      DIMENSION U(NPDE),UX(NPDE),UXX(NPDE),DFDU(NPDE,NPDE),
     *   DFDUX(NPDE,NPDE),DFDUXX(NPDE,NPDE),CMAX(1),SAVE2(1)
      COMMON /GEAR9/ EPSJ,R0,ML,MU,MW,NM1,N0ML,N0W
      ID = (IPT-1) * NPDE
      DO 40 J=1,NPDE
        UJ = U(J)
        R = EPSJ * CMAX(J)
        R = MAX(R,R0)
        U(J) = U(J) + R
        RINV = 1. / R
        CALL F(T,X,U,UX,UXX,DFDU(1,J),NPDE)
        DO 10 I=1,NPDE
   10     DFDU(I,J) = ( DFDU(I,J) - SAVE2(I+ID) ) * RINV
        U(J) = UJ
        UJ = UX(J)
        UX(J) = UX(J) + R
        CALL F(T,X,U,UX,UXX,DFDUX(1,J),NPDE)
        DO 20 I=1,NPDE
   20     DFDUX(I,J) = ( DFDUX(I,J) - SAVE2(I+ID) ) * RINV
        UX(J) = UJ
        UJ = UXX(J)
        UXX(J) = UXX(J) + R
        CALL F(T,X,U,UX,UXX,DFDUXX(1,J),NPDE)
        DO 30 I=1,NPDE
   30     DFDUXX(I,J) = ( DFDUXX(I,J) - SAVE2(I+ID) ) * RINV
        UXX(J) = UJ
   40 CONTINUE
      RETURN
      END
      SUBROUTINE EVAL(ICPT,NPDE,C,UVAL,A,ILEFT)
      IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C SUBROUTINE EVAL EVALUATES U(K), UX(K), AND UXX(K), K=1 TO NPDE,
C AT THE COLLOCATION POINT WITH INDEX ICPT USING THE VALUES OF
C THE BASIS FUNCTION COEFFICIENTS IN C AND THE BASIS FUNCTION VALUES
C STORED IN A.  THE RESULTS ARE STORED IN UVAL AS FOLLOWS..
C UVAL(K,1) = U(K), UVAL(K,2) = UX(K), AND UVAL(K,3) = UXX(K).
C
C PACKAGE ROUTINES CALLED..  NONE
C USER ROUTINES CALLED..     NONE
C CALLED BY..                GFUN,PDECOL,PSETIB
C FORTRAN FUNCTIONS USED..   NONE
C-----------------------------------------------------------------------
      DIMENSION C(NPDE,1),UVAL(NPDE,3),A(1),ILEFT(1)
      COMMON /SIZES/ NINT,KORD,IDUM(5)
      IK = ILEFT(ICPT) - KORD
      IC = 3*KORD*(ICPT-1)
      DO 10 M=1,3
        ICC = IC + KORD*(M-1)
        DO 10 J=1,NPDE
          UVAL(J,M) = 0.
          DO 10 I=1,KORD
            UVAL(J,M) = UVAL(J,M) + C(J,I+IK)*A(I+ICC)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE GFUN (T,C,UDOT,NPDE,NCPTS,A,BC,DBDU,DBDUX,DZDT,
     *                 XC,UVAL,ILEFT)
      IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C SUBROUTINE GFUN COMPUTES THE FUNCTION G(C,T), THE RIGHT-
C HAND SIDE OF THE SEMI-DISCRETE APPROXIMATION TO THE ORIGINAL
C SYSTEM OF PARTIAL DIFFERENTIAL EQUATIONS AND UPDATES THE BOUNDARY
C CONDITION INFORMATION.
C
C PACKAGE ROUTINES CALLED..  EVAL
C USER ROUTINES CALLED..     BNDRY,F
C CALLED BY..                DIFFUN,PSETIB,RES
C FORTRAN FUNCTIONS USED..   NONE
C-----------------------------------------------------------------------
      DIMENSION C(NPDE,NCPTS),UDOT(NPDE,NCPTS)
      DIMENSION A(1),BC(NPDE,NPDE,4),XC(1),UVAL(NPDE,3),ILEFT(1)
      DIMENSION DZDT(NPDE),DBDU(NPDE,NPDE),DBDUX(NPDE,NPDE)
      COMMON /SIZES/ NINT,KORD,IDUM(4),IQUAD
      DO 10 K=1,4
        DO 10 J=1,NPDE
          DO 10 I=1,NPDE
            BC(I,J,K) = 0.0
   10 CONTINUE
C-----------------------------------------------------------------------
C UPDATE THE LEFT BOUNDARY VALUES.  SAVE LEFT BOUNDARY CONDITION
C INFORMATION IN THE FIRST 2*NPDE*NPDE LOCATIONS OF BC.
C
C NOTE.. UVAL(K,1) = U(K), UVAL(K,2) = UX(K), AND UVAL(K,3) = UXX(K).
C-----------------------------------------------------------------------
      CALL EVAL(1,NPDE,C,UVAL,A,ILEFT)
      CALL BNDRY(T,XC(1),UVAL,UVAL(1,2),DBDU,DBDUX,DZDT,NPDE)
      CALL F(T,XC(1),UVAL,UVAL(1,2),UVAL(1,3),UDOT,NPDE)
      ILIM = KORD + 2
      DO 30 K=1,NPDE
        BC(K,K,1) = 1.
        IF( DBDU(K,K) .EQ. 0.0  .AND.  DBDUX(K,K) .EQ. 0.0 ) GO TO 30
        UDOT(K,1) = DZDT(K)
        DO 20 J=1,NPDE
          BC(K,J,2) = A(ILIM) * DBDUX(K,J)
          BC(K,J,1) = DBDU(K,J) - BC(K,J,2)
   20   CONTINUE
   30 CONTINUE
C-----------------------------------------------------------------------
C MAIN LOOP TO FORM RIGHT SIDE OF ODES AT THE COLLOCATION POINTS.
C-----------------------------------------------------------------------
      ILIM = NCPTS - 1
      DO 40 I=2,ILIM
        CALL EVAL(I,NPDE,C,UVAL,A,ILEFT)
        CALL F(T,XC(I),UVAL,UVAL(1,2),UVAL(1,3),UDOT(1,I),NPDE)
   40 CONTINUE
C-----------------------------------------------------------------------
C UPDATE THE RIGHT BOUNDARY VALUES.  SAVE THE RIGHT BOUNDARY CONDITION
C INFORMATION IN THE LAST 2*NPDE*NPDE LOCATIONS IN BC.
C-----------------------------------------------------------------------
      CALL EVAL(NCPTS,NPDE,C,UVAL,A,ILEFT)
      CALL F(T,XC(NCPTS),UVAL,UVAL(1,2),UVAL(1,3),UDOT(1,NCPTS),NPDE)
      CALL BNDRY(T,XC(NCPTS),UVAL,UVAL(1,2),DBDU,DBDUX,DZDT,NPDE)
      ILIM = NCPTS * 3 * KORD - KORD - 1
      DO 60 K=1,NPDE
        BC(K,K,4) = 1.
        IF( DBDU(K,K) .EQ. 0.0  .AND.  DBDUX(K,K) .EQ. 0.0 ) GO TO 60
        UDOT(K,NCPTS) = DZDT(K)
        DO 50 J=1,NPDE
          BC(K,J,3) = A(ILIM) * DBDUX(K,J)
          BC(K,J,4) = DBDU(K,J) - BC(K,J,3)
   50   CONTINUE
   60 CONTINUE
      RETURN
      END
      SUBROUTINE INTERP (TOUT, Y, N0, Y0)
      IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C SUBROUTINE INTERP COMPUTES INTERPOLATED VALUES OF THE DEPENDENT
C VARIABLE Y AND STORES THEM IN Y0.  THE INTERPOLATION IS TO THE
C POINT T = TOUT, AND USES THE NORDSIECK HISTORY ARRAY Y, AS FOLLOWS..
C                             NQ
C                  Y0(I)  =  SUM  Y(I,J+1)*S**J ,
C                            J=0
C WHERE S = -(T-TOUT)/H.
C
C PACKAGE ROUTINES CALLED..  NONE
C USER ROUTINES CALLED..     NONE
C CALLED BY..                PDECOL
C FORTRAN FUNCTIONS USED..   NONE
C-----------------------------------------------------------------------
      DIMENSION Y0(N0),Y(N0,1)
      COMMON /GEAR1/ T,H,DUMMY(4),N,IDUMMY(2),JSTART
      DO 10 I = 1,N
   10   Y0(I) = Y(I,1)
      L = JSTART + 1
      S = (TOUT - T)/H
      S1 = 1.
      DO 30 J = 2,L
        S1 = S1*S
        DO 20 I = 1,N
   20     Y0(I) = Y0(I) + S1*Y(I,J)
   30 CONTINUE
      RETURN
      END
      SUBROUTINE RES(T,H,C,V,R,NPDE,NCPTS,A,ILEFT,BC,DBDU,DBDUX,DZDT,
     *               XC,UVAL)
      IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C SUBROUTINE RES COMPUTES THE RESIDUAL VECTOR R = H*G(C,T) - A(C,T)*V
C WHERE H IS THE CURRENT TIME STEP SIZE, G IS A VECTOR, A IS A
C MATRIX, V IS A VECTOR, AND T IS THE CURRENT TIME.
C
C PACKAGE ROUTINES CALLED..  GFUN
C USER ROUTINES CALLED..     NONE
C CALLED BY..                STIFIB
C FORTRAN FUNCTIONS USED..   NONE
C-----------------------------------------------------------------------
      DIMENSION C(NPDE,NCPTS),R(NPDE,1),V(NPDE,1)
      DIMENSION A(1),ILEFT(1),BC(NPDE,NPDE,4),XC(1),UVAL(1)
      DIMENSION DBDU(NPDE,NPDE),DBDUX(NPDE,NPDE),DZDT(NPDE)
      COMMON /SIZES/ NINT,KORD,NCC,IDUM(3),IQUAD
C-----------------------------------------------------------------------
C FORM G(C,T) AND STORE IN R.
C-----------------------------------------------------------------------
      CALL GFUN(T,C,R,NPDE,NCPTS,A,BC,DBDU,DBDUX,DZDT,XC,UVAL,ILEFT)
C-----------------------------------------------------------------------
C FORM THE FIRST AND LAST BLOCK ROWS OF THE RESIDUAL VECTOR
C WHICH ARE DEPENDENT ON THE BOUNDARY CONDITIONS.
C-----------------------------------------------------------------------
      ILIM = NCPTS - 1
      DO 20 I=1,NPDE
        SUM1 = 0.0
        SUM2 = 0.0
        DO 10 J=1,NPDE
          SUM1 = SUM1 + BC(I,J,1) * V(J,1) + BC(I,J,2) * V(J,2)
          SUM2 = SUM2 + BC(I,J,3) * V(J,ILIM) + BC(I,J,4) * V(J,NCPTS)
   10   CONTINUE
        R(I,1) = H * R(I,1) - SUM1
        R(I,NCPTS) = H * R(I,NCPTS) - SUM2
   20 CONTINUE
C-----------------------------------------------------------------------
C FORM THE REMAINING COMPONENTS OF THE RESIDUAL VECTOR.
C-----------------------------------------------------------------------
      DO 50 ICPTS=2,ILIM
        I2 = (ICPTS-1) * KORD * 3
        ICOL = ILEFT(ICPTS) - KORD
        DO 40 JJ=1,NPDE
          SUM1 = 0.
          DO 30 J=1,KORD
            SUM1 = SUM1 + A(I2+J) * V(JJ,ICOL+J)
   30       CONTINUE
          R(JJ,ICPTS) = H*R(JJ,ICPTS) - SUM1
   40     CONTINUE
   50   CONTINUE
      RETURN
      END
      SUBROUTINE VALUES(X,USOL,SCTCH,NDIM1,NDIM2,NPTS,NDERV,WORK)
      IMPLICIT REAL*8(A-H,O-Z)
C-----------------------------------------------------------------------
C SUBROUTINE VALUES COMPUTES THE SOLUTION U AND THE FIRST NDERV
C DERIVATIVES OF U AT THE NPTS POINTS X AND AT TIME TOUT AND RETURNS
C THEM IN THE ARRAY USOL.  THIS ROUTINE MUST BE USED TO OBTAIN
C SOLUTION VALUES SINCE PDECOL DOES NOT RETURN ANY SOLUTION VALUES
C TO THE USER.  SEE PDECOL.
C
C THE CALLING PARAMETERS ARE...
C   X     =  AN ARBITRARY VECTOR OF SPATIAL POINTS OF LENGTH NPTS AT
C            WHICH THE SOLUTION AND THE FIRST NDERV DERIVATIVE VALUES
C            ARE TO BE CALCULATED.  IF X .LT. XLEFT ( X .GT. XRIGHT )
C            THEN THE PIECEWISE POLYNOMIAL OVER THE LEFTMOST ( RIGHT-
C            MOST ) INTERVAL IS EVALUATED TO CALCULATE THE SOLUTION
C            VALUES AT THIS UNUSUAL VALUE OF  X.  SEE PDECOL.
C
C   USOL  =  AN ARRAY WHICH CONTAINS THE SOLUTION AND THE FIRST
C            NDERV DERIVATIVES OF THE SOLUTION AT ALL THE POINTS IN
C            THE INPUT VECTOR X.  IN PARTICULAR, USOL(J,I,K) CONTAINS
C            THE VALUE OF THE (K-1)-ST DERIVATIVE OF THE J-TH PDE
C            COMPONENT AT THE I-TH POINT OF THE X VECTOR FOR
C            J = 1 TO NPDE, I = 1 TO NPTS, AND K = 1 TO NDERV+1.
C
C   SCTCH =  A USER SUPPLIED WORKING STORAGE ARRAY OF LENGTH AT LEAST
C            KORD*(NDERV+1).  SEE BELOW AND PDECOL FOR DEFINITIONS OF
C            THESE PARAMETERS.
C
C   NDIM1 =  THE FIRST DIMENSION OF THE OUTPUT ARRAY USOL IN THE CALLING
C            PROGRAM.  NDIM1 MUST BE .GE. NPDE.
C
C   NDIM2 =  THE SECOND DIMENSION OF THE OUTPUT ARRAY USOL IN THE
C            CALLING PROGRAM.  NDIM2 MUST BE .GE. NPTS.
C
C   NPTS  =  THE NUMBER OF POINTS IN THE X VECTOR.
C
C   NDERV =  THE NUMBER OF DERIVATIVE VALUES OF THE SOLUTION THAT ARE
C            TO BE CALCULATED.  NDERV SHOULD BE LESS THAN KORD SINCE
C            THE KORD-TH DERIVATIVE OF A POLYNOMIAL OF DEGREE KORD-1
C            IS EQUAL TO ZERO.  SEE PDECOL.
C
C   WORK  =  THE WORKING STORAGE ARRAY DESCRIBED IN PDECOL.
C
C PACKAGE ROUTINES CALLED..  BSPLVD,INTERV
C USER ROUTINES CALLED..     NONE
C CALLED BY..                USERS MAIN PROGRAM
C FORTRAN FUNCTIONS USED..   NONE
C
C-----------------------------------------------------------------------
      DIMENSION USOL(NDIM1,NDIM2,NDERV),X(NPTS),SCTCH(1),WORK(1)
      COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD
      COMMON /ISTART/ IW1,IW2,IW3,IW4,IW5,IW6,IDUM(12)
      DATA ILEFT/0/, MFLAG/0/
      NDERV1 = NDERV + 1
      DO 20 IPTS=1,NPTS
        CALL INTERV(WORK(IW2),NCPTS,X(IPTS),ILEFT,MFLAG)
        CALL BSPLVD(WORK(IW2),KORD,X(IPTS),ILEFT,SCTCH,NDERV1)
        IK = ILEFT - KORD
        DO 10 M=1,NDERV1
          I1 = (M-1)*KORD
          DO 10 K=1,NPDE
            USOL(K,IPTS,M) = 0.
            DO 10 I=1,KORD
              I2 = (I+IK-1)*NPDE + IW6 - 1
              USOL(K,IPTS,M) = USOL(K,IPTS,M) + WORK(I2+K) * SCTCH(I+I1)
   10   CONTINUE
   20 CONTINUE
      RETURN
      END
