! Copyright (c) 2006, Paul Muir and Larry Shampine.
! Paul Muir, Mathematics and Computing Science, Saint Mary's University.
! Larry Shampine, Mathematics, Southern Methodist University.
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!     * Redistributions of source code must retain the above copyright
!       notice, this list of conditions and the following disclaimer.
!     * Redistributions in binary form must reproduce the above copyright
!       notice, this list of conditions and the following disclaimer in the
!       documentation and/or other materials provided with the distribution.
!     * Neither Saint Mary's University nor Southern Methodist University nor 
!       the names of its contributors may be used to endorse or promote products
!       derived from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY Paul Muir and Larry Shampine ''AS IS'' AND ANY
! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
! WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL Paul Muir and Larry Shampine BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
! ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

! L.F. Shampine, P.H. Muir, H. Xu, A user-friendly Fortran BVP solver, 
! J. Numer. Anal. Indust. Appl. Math., 1, 2006, 201--217.

  MODULE BVP_M
! In the declarations part of this module, we first define the BVP_SOL type, 
! which represents the numerical solution. We next define the interface 
! BVP_INIT, which is overloaded to provide a number of functions for allowing 
! the user to specify initial information for the computation, the interface 
! BVP_EXTEND, which is overloaded to provide functions for extending the 
! numerical solution to a new domain that are useful in the context of parameter
! continuation, and the interface BVP_EVAL, which is overloaded to provide a
! function for retrieving solution information. Next a number of global variables 
! employed throughout the code are defined.
!
! A number of the PUBLIC routines available to the user access the BVP_SOL 
! structured type and allocate or deallocate the array fields of this type. 
! In order to avoid "memory leaks", care must be taken to explicitly
! deallocate these array fields prior to reallocating them and at the end of
! the computation. This is particularly important when BVP_SOLVER is called 
! several times within the same program. The following table provides pre and 
! post conditions on array fields of SOL for public routines that access SOL. 
! For each routine the pre condition specifies which array fields of SOL must 
! be allocated prior to calling the routine, while the post condition specifies 
! which array fields of SOL will be allocated upon return from the routine.
! 
! BVP_SOLVER     Pre: X, Y,(optionally) PARAMETERS.
!                Post/Successful Return: X, Y,(optionally) PARAMETERS,IWORK,WORK.
!                Post/Unsuccessful Return: No fields of SOL are allocated.
!
! BVP_EVAL       Pre: X, Y,(optionally) PARAMETERS, IWORK, WORK.
!                Post: X, Y,(optionally) PARAMETERS, IWORK, WORK.
!  
! BVP_INIT       Pre: No fields of SOL are allocated.
!                Post: X, Y,(optionally) PARAMETERS.
!
! BVP_EXTEND     Pre: X, Y,(optionally) PARAMETERS, IWORK, WORK.
!                Post: X, Y,(optionally) PARAMETERS.
!
! BVP_SAVE       Pre: X, Y,(optionally) PARAMETERS, IWORK, WORK.
!                Post: X, Y,(optionally) PARAMETERS, IWORK, WORK.
!
! BVP_GET        Pre: No fields of SOL are allocated.
!                Post: X, Y,(optionally) PARAMETERS, IWORK, WORK.
!
! BVP_TERMINATE  Pre: X, Y,(optionally) PARAMETERS, IWORK, WORK.
!                Post: No fields of SOL are allocated.
!-------------------------------------------------------------------------------

  IMPLICIT NONE
  
! Declare everything in the module to be PRIVATE unless specifically declared
! to be PUBLIC.
  PRIVATE
  
! PUBLIC subroutines/functions defined in this module.  
  PUBLIC :: BVP_SOLVER,BVP_EVAL,BVP_INIT,BVP_EXTEND,BVP_LINSPACE,            &
            BVP_SAVE,BVP_GET,BVP_TERMINATE

!----------------------------------------------------------------------------
! Type definitions:
! Set up structured type that will store solution and associated information.
  TYPE, PUBLIC :: BVP_SOL

    INTEGER :: NODE,NPAR ! number of ODEs; number of unknown parameters.
    INTEGER :: LEFTBC ! number of left boundary conditions.
    INTEGER :: NPTS,INFO ! number of points in current mesh; indicator of
                         ! success (INFO=0) or failure (INFO=-1) of computation.
    INTEGER :: MXNSUB ! maximum number of subintervals allowed.

    DOUBLE PRECISION, DIMENSION(:), POINTER :: X ! current mesh.

    DOUBLE PRECISION, DIMENSION(:,:), POINTER :: Y ! current solution - Ith 
!   column is solution approximation at Ith mesh point in first NODE locations; 
!   last NPAR locations contain approximation to unknown parameters.

    DOUBLE PRECISION, DIMENSION(:), POINTER :: PARAMETERS ! Unknown parameters.

    INTEGER, DIMENSION(:), POINTER :: IWORK ! integer workspace.
    DOUBLE PRECISION, DIMENSION(:), POINTER :: WORK ! double precision workspace.
!   The workspaces are used to save information that can be accessed by SOL_EVAL 
!   for evaluation of the numerical solution.

  END TYPE BVP_SOL
  
!-------------------------------------------------------------------------------
! Interface definitions:

  INTERFACE BVP_INIT
! The purpose of each of the functions available in this interface is to set 
! initial values for the fields of the solution structure SOL and optionally 
! to set the global variable, MXNSUB, the maximum number of subintervals, based 
! on the value of the MAX_NUM_SUBINTERVALS argument. This interface consists of
! three functions, GUESS_1, GUESS_2, and GUESS_3, which give the user with three 
! options for specifying the initial guess for the solution. Through each function 
! the user also provides values for NODES, the number of ODEs, LEFTBC, the number 
! of left boundary conditions, and X, the initial mesh. (The SIZE of X is used to 
! set the NPTS field, the number of mesh points.) Optionally an initial guess for 
! unknown parameters can be provided through P. The INFO field of SOL is set to 0. 
! If X contains only end points, a default mesh of 10 equally spaced points is set.

    MODULE PROCEDURE GUESS_1,GUESS_2,GUESS_3
    
!   GUESS_?(NODE,LEFTBC,X,?,P,MAX_NUM_SUBINTERVALS) RESULT(SOL)

!   GUESS_1 - the ? argument is a vector of NODE components which represents the 
!            constant vector to be assigned to all columns of the Y field of SOL,
!   GUESS_2 - the ? argument is an array of size NODE by NPTS representing, in 
!            each column, an initial solution guess at the corresponding mesh 
!            point,
!   GUESS_3 - the ? argument is the name of a subroutine that can be called to 
!            give a solution guess vector of NODE components when it is called
!            with a scalar value X. This function has the form FCN(X,Y). GUESS_3
!            calls FCN to get solution guesses for each of the points in the
!            array X and loads them in corresponding columns of the Y field.

  END INTERFACE ! BVP_INIT


  INTERFACE BVP_EXTEND
! The purpose of each of the functions available within this interface
! is to extend a solution structure on [A,B] to a larger interval [ANEW,BNEW].

    MODULE PROCEDURE EXTEND_1,EXTEND_2

!   FUNCTION EXTEND_1(SOLIN,ANEW,BNEW,YANEW,YBNEW,P,MAX_NUM_SUBINTERVALS) &
!                                                         RESULT(SOLOUT)   
!   This function allows the solution in a structure of type BVP_SOL to be
!   extended to a solution on a bigger interval. A typical call has the form

!       SOLOUT = BVP_EXTEND(SOLIN,ANEW,BNEW,YANEW,YBNEW)

!   Here SOLIN is the input solution structure and SOLOUT is the output 
!   structure.  If the interval used to compute SOLIN was [A,B], the new 
!   guess is for [ANEW,BNEW]. If ANEW < A, the mesh is extended to ANEW and 
!   Y is extended to the value YANEW there. If ANEW >= A, the arguments ANEW
!   and YANEW are ignored. The right end of the interval is treated in the 
!   same way. If present, SOLIN%PARAMETERS is used in SOLOUT. To supply a 
!   different guess P for unknown parameters, use the optional argument P.  
!   To increase the global variable, MXNSUB, the maximum number of subintervals 
!   allowed, use the optional argument MAX_NUM_SUBINTERVALS.


!   FUNCTION EXTEND_2(SOLIN,ANEW,BNEW,ORDER,P,MAX_NUM_SUBINTERVALS) &
!                                                          RESULT(SOLOUT)   
!   This function allows the solution in a structure of type BVP_SOL to be
!   extended to a solution on a bigger interval. A typical call has the form

!       SOLOUT = BVP_EXTEND(SOLIN,ANEW,BNEW,ORDER)

!   Here SOLIN is the input solution structure and SOLOUT is the output 
!   structure.  If the interval used to compute SOLIN was [A,B], the new 
!   guess is for [ANEW,BNEW]. If ANEW < A, the mesh is extended to ANEW and 
!   Y is extended to a value YANEW there. If ANEW >= A, no action is taken
!   at the left end.  The value YANEW is computed by extrapolation.  The 
!   optional integer ORDER indicates the order of the polynomial used.  
!   If ORDER = 1, linear extrapolation is done and otherwise, constant 
!   extrapolation (ORDER = 0). The right end of the interval is treated in 
!   the same way. If present, SOLIN%PARAMETERS is used in SOLOUT. To supply 
!   a different guess P for unknown parameters, use the optional argument P.  
!   To increase the global variable, MXNSUB, the maximum number of subintervals 
!   allowed, use the optional argument MAX_NUM_SUBINTERVALS.

  END INTERFACE ! BVP_EXTEND


  INTERFACE BVP_EVAL
! Assuming a successful computation, this interface provides three routines
! for evaluation of solution information. EVAL_S evaluates the numerical
! solution at a single point. EVAL_V evaluates the numerical solution at a 
! set of points. EVAL_P returns values for unknown parameters, if there are any.
    
    MODULE PROCEDURE EVAL_S,EVAL_V,EVAL_P

!   SUBROUTINE EVAL_S(SOL,T,Z,DERIVATIVE) 
!   Evaluate solution (and optionaly first derivative) at a single point
!   T. The solution is returned in Z; the derivative in DERIVATIVE if 
!   this parameter is present. 

!   SUBROUTINE EVAL_V(SOL,T,Z,DERIVATIVE)
!   Evaluate solution (and optionally first derivative) at all the points
!   of the vector T. Solutions are returned in the columns of the array Z; 
!   derivatives are returned in the columns of the array DERIVATIVE if
!   this parameter is present. 

!   SUBROUTINE EVAL_P(SOL,PARAMETERS)
!   Return parameter values in PARAMETERS. 

  END INTERFACE ! BVP_EVAL

!-------------------------------------------------------------------------------
! Global Constants and Variables:

  DOUBLE PRECISION, PARAMETER :: UROUND=EPSILON(1D0) ! Unit roundoff.
  
  DOUBLE PRECISION, PARAMETER :: LAMBDA_MIN=0.01D0 ! Minimum value for damping 
                                                   ! factor for damped Newton 
                                                   ! iteration.
  
  DOUBLE PRECISION :: SQRTU ! Square root of unit roundoff.

  DOUBLE PRECISION, DIMENSION(:,:), POINTER :: MSING,FSING,BCSING
  LOGICAL :: SINGULAR  ! Indicate whether solving a singular BVP.


! Declaration of constants, variables and arrays associated with discrete
! and continuous Runge-Kutta formulas.

  INTEGER, PARAMETER :: MXS=10 ! Maximum number of stages of Runge-Kutta method.

  INTEGER :: C_S, C_S_STAR, C_P
! C_S, number of discrete stages; C_S_STAR, total number of stages required to 
! form the interpolant on each subinterval. It includes all the stages of the 
! discrete formula plus the additional stages required for the interpolant; 
! C_P is the order of the interpolant => local error is O(h**(p+1)), where h 
! is the subinterval size.

  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: C_C, C_V, C_B
  DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: C_X
! C_C, C_V, C_B, and C_X are coefficients that define the discrete Runge-Kutta 
! formula.  C_C, C_V, and C_B correspond, respectively, to the Runge-Kutta 
! coefficient vectors c, v, and b; C_X corresponds to the Runge-Kutta coefficient 
! matrix X.

  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: C_C_STAR, C_V_STAR
  DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: C_X_STAR
! C_C_STAR, C_V_STAR, and C_X_STAR are coefficients that define the continuous 
! Runge-Kutta formula.  C_C_STAR and C_V_STAR correspond, respectively, to the 
! Runge-Kutta coefficient vectors c^star and v^star; C_X_STAR corresponds to 
! the Runge-Kutta coefficient matrix X^star.

  DOUBLE PRECISION :: C_TAU_STAR ! Relative position of one of the sample points 
                                 ! within each subinterval. The other sample 
                                 ! point is 1-C_TAU_STAR.


  INTEGER :: MXNSUB ! maximum number of subintervals --  this is used to provide 
! an upper bound that can be checked in order to determine when the next mesh 
! would too large - this is is preferable to termination based on a failed memory 
! allocation call. (Recommended by Skip Thompson, May 26th, 2005.) The default 
! value is 3000, but it can be changed with the option MAX_NUM_SUBINTERVALS in 
! BVP_INIT and BVP_EXTEND.

  INTEGER :: PROFILE ! Global variable to specify level of TRACE. 

  LOGICAL :: G_STOP_ON_FAIL ! Global variable to allow access to STOP_ON_FAIL
                            ! argument to BVP_SOLVER.

  DOUBLE PRECISION :: BVP_TOL, NEWTON_TOL ! User tolerance for size of defect of 
! approximate solution and tolerance for Newton iteration.  BVP_TOL is input TOL 
! or a default value.

  INTEGER :: NPTS,NSUB  ! Number of points and number of subintervals for current
                        ! mesh.

  LOGICAL :: HAVE_DFDY,HAVE_DBCDY ! Record presence of DFDY and DBCDY functions.

  INTEGER :: NODE,NPAR,NEQN,NEQNSQ,BVP_METHOD,LEFTBC,RIGHTBC
! NODE, number of ODEs, NPAR, number of unknown parameters; number of equations 
! is NEQN = NODE+NPAR. NEQNSQ = NEQN**2. BVP_METHOD identifies the Runge-Kutta 
! method to be used and LEFTBC and RIGHTBC are the number of left and right BCs, 
! respectively. 

  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: G_X, G_Y
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: G_K_DISCRETE
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: G_K_INTERP
! Allocatable arrays to be used for storage of the mesh, discrete solution, and 
! discrete and continuous Runge-Kutta stages, respectively.

CONTAINS

  FUNCTION BVP_SOLVER(SOL_INIT,FSUB,BCSUB,SINGULARTERM,METHOD,TOL,DFDY,          &
                      DBCDY,TRACE,STOP_ON_FAIL,YERROR) RESULT(SOL)
! This function is called by the user's main program to solve the boundary 
! value problem. It assumes that BVP_INIT has been called previously to form
! the guess structure, SOL_INIT. The ODEs and boundary conditions are made 
! available through the user-provided routines, FSUB and BCSUB, and the
! corresponding partial derivatives are optionally available through the 
! user-provided routines, DFDY and DBCDY. If there is a singular term, it is
! specified by the matrix SINGULARTERM. The Runge-Kutta methods to be employed 
! are specified by the value of METHOD.  The possibilities are 2 for the second 
! order method, 4 for the fourth order method, 6 for the sixth order method. 
! The default is 4. The user-requested tolerance is provided through TOL. The
! default is 1D-6.  The level of monitoring of the computation is specified by 
! TRACE. The default is 0 for no output unless the computation fails, 1 for 
! intermediate output, 2 for full output; for these values of TRACE the solver
! STOPs when there is a failure. If TRACE = -1, there is no output under any 
! circumstances and the solver does not STOP when there is a failure--the user 
! must check SOL%INFO to determine if computation was successful. The user can 
! set STOP_ON_FAIL to .FALSE. if the code should not STOP on failure; the default 
! (STOP_ON_FAIL=.TRUE.) is to STOP if the computation fails. This code controls
! the error in the ODEs and the boundary conditions.  Providing the optional
! argument YERROR causes the code to estimate the error in the solution itself.
! If present, the code will estimate the maximum error in a solution component 
! y(i) relative to 1 + |y(i)| for all i and all mesh points and return this
! quantity as YERROR. The same estimate applies to unknown parameters.
!-------------------------------------------------------------------------------
!     CALLED BY: User program
!     CALLS: RK_TABLEAU,INTERP_TABLEAU,CHECK_STAT,MIRKDC,EYE,DUMMY_DF,DUMMY_DBC,
!            PINVSQ
!-------------------------------------------------------------------------------
!   Input arguments:
    TYPE(BVP_SOL) :: SOL_INIT
    DOUBLE PRECISION, DIMENSION(:,:), OPTIONAL :: SINGULARTERM
    INTEGER, OPTIONAL :: METHOD,TRACE
    LOGICAL, OPTIONAL :: STOP_ON_FAIL
    DOUBLE PRECISION, OPTIONAL :: TOL
    
!   User-supplied subroutines:
    EXTERNAL FSUB,BCSUB,DFDY,DBCDY
    OPTIONAL :: DFDY,DBCDY
           
!   Output arguments:
    TYPE(BVP_SOL) ::  SOL
    DOUBLE PRECISION, OPTIONAL :: YERROR
       
!   Local variables:
    INTEGER :: N_IWORK,N_WORK ! Dimensions of work array fields of SOL, IWORK 
                              ! and WORK.
    INTEGER :: IER ! Dymanic array allocation error flag.
    INTEGER :: NDEX,M  ! Array index and loop index, respectively.
    INTEGER :: INFO ! Communication on status of computation (error flag) from
                    ! MIRKDC: INFO=0 for success; INFO=-1 for failure.
    INTEGER :: I, GE_INFO ! Loop index; INFO flag for global error NewIter call.
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: TEMPY  ! Temporary storage.
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: ORIGY,WT ! Temporary storage.

!   Global variables:
!   NODE is the number of ODEs.
!   NPAR is the number of unknown parameters.
!   NEQN equals NODE + NPAR.
!   NEQNSQ equals NEQN**2.
!   LEFTBC and RIGHTBC are the number of left and right boundary conditions, 
!     respectively. 
!   MXNSUB is the maximum number of subintervals.
!   MSING,FSING,BCSING, and SINGULAR are quantities used in dealing with 
!     singular BVPs.
!   UROUND is unit roundoff.
!   SQRTU is the square root of unit roundoff.
!   BVP_METHOD identifies which Runge-Kutta methods are to be employed. 
!   C_S and C_C_STAR are the number of discrete Runge-Kutta stages and continuous 
!     Runge-Kutta stages, respectively.     
!   C_C, C_V, C_B, and C_X are coefficients that define the discrete Runge-Kutta 
!     formula.  
!   C_C_STAR, C_V_STAR, and C_X_STAR are coefficients that define the continuous 
!     Runge-Kutta formula.  
!   BVP_TOL is input TOL or the default value.
!   NEWTON_TOL is the tolerance applied to the Newton iteration for the 
!     determination of the discrete solution.
!   PROFILE gives the level of output to be provided.
!   NPTS is the number of points for the current mesh.
!   NSUB is the number of subintervals for the current mesh.
!   G_X, current mesh; its dimension is 0:NSUB.
!   G_Y, current discrete solution approximation, in column form; 
!     its dimension is NEQN*(NSUB+1).
!   G_K_DISCRETE contains the discrete Runge-Kutta stages.  It has dimension
!     C_S*NEQN*NSUB.
!   G_K_INTERP contains the discrete continuous Runge-Kutta stages. It has 
!     dimension (C_S_STAR-C_S)*NEQN*NSUB.
!   HAVE_DFDY,HAVE_DBCDY record presence of DFDY and DBCDY functions.
!   G_STOP_ON_FAIL determines whether or not the code stops itself after a 
!     failure or simply returns with INFO=-1.
!------------------------------------------------------------------------

!   Load values for the global variables, NODE, NPAR, and LEFTBC from SOL_INIT.
!   Compute values for global variables NEQN, NEQNSQ, and RIGHTBC.
    NODE = SOL_INIT%NODE
    NPAR = SOL_INIT%NPAR
    NEQN = NODE + NPAR
    NEQNSQ = NEQN**2
    LEFTBC = SOL_INIT%LEFTBC
    RIGHTBC = NEQN - LEFTBC

!   Load value for global variable MXNSUB from SOL_INIT.
    MXNSUB = SOL_INIT%MXNSUB    

!   Set global variable INFO to initial value defined in SOL_INIT.
    INFO = SOL_INIT%INFO

!   Set global variables NPTS, NSUB.
    NPTS = SOL_INIT%NPTS
    NSUB = NPTS - 1 

!   If the problem has a singular term, set up associated quantities that 
!   will be needed for the computation.
    IF (PRESENT(SINGULARTERM)) THEN
      SINGULAR = .TRUE.                 
      ALLOCATE(MSING(NODE,NODE),FSING(NODE,NODE),BCSING(NODE,NODE),STAT=IER)
      CALL CHECK_STAT(IER)
      MSING = SINGULARTERM
      FSING = PINVSQ( EYE(NODE)-MSING ,NODE)
      BCSING = EYE(NODE) - MATMUL(PINVSQ(MSING,NODE), MSING)
    ELSE
      SINGULAR = .FALSE.
    END IF

!   Set global variable SQRTU.
    SQRTU = SQRT(UROUND)
      
!   METHOD is optional and the default value is 4. METHOD is communicated 
!   throughout this code by the global variable BVP_METHOD. 
    IF (PRESENT(METHOD)) THEN
      BVP_METHOD = METHOD
    ELSE
      BVP_METHOD = 4
    END IF
      
!   Define C_S and C_S_STAR values based on BVP_METHOD.
    SELECT CASE (BVP_METHOD)
      CASE (2)
        C_S = 2
        C_S_STAR = 2
      CASE (4)
        C_S = 3
        C_S_STAR = 4
      CASE (6)
        C_S = 5
        C_S_STAR = 8
    END SELECT 

!   Define remaining discrete Runge-Kutta coefficients.
    CALL RK_TABLEAU
!   Define remaining continuous Runge-Kutta coefficients.
    CALL INTERP_TABLEAU

!   Set global variables BVP_TOL and NEWTON_TOL.
    IF (PRESENT(TOL)) THEN
!     IF (TOL < 100*UROUND) THEN
!       PRINT *,' Input tolerance is too small; choose a larger value.'
!       STOP
!     END IF       
      BVP_TOL = TOL
    ELSE
      BVP_TOL = 1D-6
    END IF
    NEWTON_TOL = BVP_TOL/100D0
    
!   Set level of output with global variable PROFILE, based on value of TRACE.
    IF (PRESENT(TRACE)) THEN
      IF (TRACE < -1 .OR. TRACE > 2) THEN
        PRINT *,' TRACE must be -1, 0, 1, or 2.'
        STOP
      END IF
      PROFILE = TRACE
    ELSE
      PROFILE = 0
    ENDIF
     
!   Set global variable G_STOP_ON_FAIL. Default is to stop on a failure.
    IF (PRESENT(STOP_ON_FAIL)) THEN
      G_STOP_ON_FAIL = STOP_ON_FAIL
    ELSE
      G_STOP_ON_FAIL = .TRUE.
    ENDIF

!   Allocate space for global arrays G_X, G_Y, G_K_DISCRETE and G_K_INTERP; 
!   these will be used within MIRKDC, where they will be deallocated and 
!   reallocated, as the mesh is refined and the solution approximation improved. 
    ALLOCATE(G_X(0:NSUB),G_Y(NPTS*NEQN),G_K_DISCRETE(C_S*NEQN*NSUB),             &
             G_K_INTERP((C_S_STAR-C_S)*NEQN*NSUB),STAT=IER)
    CALL CHECK_STAT(IER) 

!   Store the initial mesh - obtained from SOL_INIT - in G_X.
    G_X = SOL_INIT%X

!   Set up a copy of initial solution in G_Y. Account for unknown parameters.
    IF (NPAR == 0) THEN
      G_Y = RESHAPE(SOL_INIT%Y,(/NPTS*NEQN/))
    ELSE
      ALLOCATE(TEMPY(NEQN,NPTS),STAT=IER)
      CALL CHECK_STAT(IER)
      TEMPY(1:NODE,1:NPTS) = SOL_INIT%Y
      DO M = 1,NPTS
        TEMPY(NODE+1:NEQN,M) = SOL_INIT%PARAMETERS
      END DO
      G_Y = RESHAPE(TEMPY,(/NPTS*NEQN/))
      DEALLOCATE(TEMPY,STAT=IER)
      CALL CHECK_STAT(IER)
    END IF
     
!   If problem includes a singular term, impose necessary condition on  
!   approximation at left end point.
    IF (SINGULAR) THEN
      G_Y(1:NODE) = MATMUL(BCSING,G_Y(1:NODE))
    END IF

!   All information has now been copied from SOL_INIT. Deallocate array 
!   fields in SOL_INIT: X, Y, and optionally PARAMETERS.

    DEALLOCATE(SOL_INIT%X, SOL_INIT%Y,STAT=IER)
    CALL CHECK_STAT(IER)
    IF (NPAR > 0) THEN
      DEALLOCATE(SOL_INIT%PARAMETERS,STAT=IER)
      CALL CHECK_STAT(IER)
    END IF
    
!   Call MIRKDC with argument list that depends on whether or not 
!   DFDY and DGDY are present.    
    HAVE_DFDY = PRESENT(DFDY)
    HAVE_DBCDY = PRESENT(DBCDY)
    IF (HAVE_DFDY) THEN
      IF (HAVE_DBCDY) THEN
        CALL MIRKDC(INFO,FSUB,BCSUB,DFDY,DBCDY)
      ELSE
        CALL MIRKDC(INFO,FSUB,BCSUB,DFDY,DUMMY_DBC)
      END IF  
    ELSE
      IF (HAVE_DBCDY) THEN
        CALL MIRKDC(INFO,FSUB,BCSUB,DUMMY_DF,DBCDY)
      ELSE
        CALL MIRKDC(INFO,FSUB,BCSUB,DUMMY_DF,DUMMY_DBC)
      END IF
    END IF
     
!   Failure of the MIRKDC call is indicated by INFO = -1
!   (Maximum mesh size, MXNSUB, exceeded during attempt to
!   determine new mesh.)
!   Load INFO and NPTS fields of SOL (even on failure.)
    SOL%INFO = INFO
    NPTS = NSUB + 1
    SOL%NPTS = NPTS

!   If computation performed by MIRKDC was successful ...
    IF (INFO >= 0) THEN

!     Copy NODE, NPAR, LEFTBC, MXNSUB values from global 
!     variables to SOL fields.
      SOL%NODE = NODE
      SOL%NPAR = NPAR
      SOL%LEFTBC = LEFTBC
      SOL%MXNSUB = MXNSUB 
      
!     Load mesh and discrete solution into SOL.

      ALLOCATE(SOL%X(NPTS),SOL%Y(NEQN,NPTS),STAT=IER)
      CALL CHECK_STAT(IER)
      SOL%X = G_X
      SOL%Y = RESHAPE(G_Y,(/NEQN,NPTS/))
     
!     If present, load parameters into SOL.
      IF (NPAR > 0) THEN
        ALLOCATE(SOL%PARAMETERS(NPAR),STAT=IER)
        CALL CHECK_STAT(IER)
        SOL%PARAMETERS = SOL%Y(NODE+1:NEQN,1)
        ALLOCATE(TEMPY(NODE,NPTS),STAT=IER)
        CALL CHECK_STAT(IER)
        TEMPY = SOL%Y(1:NODE,1:NPTS)
        DEALLOCATE(SOL%Y)
        CALL CHECK_STAT(IER)
        ALLOCATE(SOL%Y(NODE,NPTS),STAT=IER)
        CALL CHECK_STAT(IER)
        SOL%Y = TEMPY
        DEALLOCATE(TEMPY)
        CALL CHECK_STAT(IER)
      END IF

!     Load SOL%IWORK and SOL%WORK fields needed by BVP_EVAL. 
      N_IWORK = 3
      N_WORK = C_S_STAR*NEQN*NSUB + (NSUB+1) + (NSUB+1)*NEQN

      ALLOCATE(SOL%IWORK(N_IWORK),SOL%WORK(N_WORK),STAT=IER)
      CALL CHECK_STAT(IER)
 
!     Save global variable values needed by BVP_EVAL. 
      SOL%IWORK(1) = C_S
      SOL%IWORK(2) = C_S_STAR
      SOL%IWORK(3) = BVP_METHOD

!     Copy G_K_DISCRETE and G_K_INTERP into first two sections of WORK array.
      SOL%WORK(1:C_S*NEQN*NSUB) = G_K_DISCRETE 
      NDEX = C_S*NEQN*NSUB
      SOL%WORK(NDEX+1:NDEX+(C_S_STAR-C_S)*NEQN*NSUB) = G_K_INTERP

!     Copy G_X into the third part of SOL%WORK           
      NDEX = C_S_STAR*NEQN*NSUB
      SOL%WORK(NDEX+1:NDEX+NSUB+1) = G_X

!     Copy G_Y into the fourth part of SOL%WORK
      NDEX = C_S_STAR*NEQN*NSUB + (NSUB+1)
      SOL%WORK(NDEX+1:NDEX+(NSUB+1)*NEQN) = G_Y

      IF (PRESENT(YERROR)) THEN
!       Estimate the error in solution component y(i) relative to 1 + |y(i)|.
!       Return the maximum over all components and mesh points as YERROR. The
!       same estimate applies to unknown parameters.  The estimate is computed 
!       by extrapolation from the current solution SOL to a solution on the 
!       "doubled" mesh.

!       We will need storage for the doubled mesh, the corresponding solution,
!       and the corresponding discrete Runge-Kutta stages. Deallocate current 
!       working storage arrays and reallocate them for the doubled mesh.
        DEALLOCATE(G_X,G_Y,G_K_DISCRETE,STAT=IER)
        CALL CHECK_STAT(IER)
        ALLOCATE(G_X(0:2*NSUB),G_Y((2*NSUB+1)*NEQN),                             &
                 G_K_DISCRETE(C_S*NEQN*(2*NSUB)),STAT=IER)
        CALL CHECK_STAT(IER)
        
!       Initialize new mesh, G_X, and set initial guess, G_Y, based
!       on original solution.
        G_X(0:2*NSUB:2) = SOL%X
        G_X(1:2*NSUB-1:2) = 0.5D0*(G_X(0:2*NSUB-2:2) + G_X(2:2*NSUB:2))
      
        DO I = 1,2*NSUB+1
!         Solution values.
          CALL BVP_EVAL(SOL,G_X(I-1),G_Y((I-1)*NEQN+1:(I-1)*NEQN+NODE))
!         Parameter values. 
          IF (NPAR > 0 ) THEN
            CALL BVP_EVAL(SOL,G_Y((I-1)*NEQN+NODE+1:(I-1)*NEQN+NEQN))
          END IF
        END DO
        
!       Reset global variables NSUB and NPTS for doubled mesh.
        NSUB = 2*NSUB
        NPTS = NSUB+1
        
!       Reset NEWTOL to be slightly sharper for a more accurate computation.
        NEWTON_TOL = NEWTON_TOL/2D0
        
!       Setup and solve the Newton system asscoiated with the doubled mesh.
        IF (HAVE_DFDY) THEN
          IF (HAVE_DBCDY) THEN
            CALL NEWITER(NEQN,NSUB,G_X,G_Y,GE_INFO,G_K_DISCRETE,             &
                         FSUB,BCSUB,DFDY,DBCDY)
          ELSE
            CALL NEWITER(NEQN,NSUB,G_X,G_Y,GE_INFO,G_K_DISCRETE,             &
                         FSUB,BCSUB,DFDY,DUMMY_DBC)
          END IF  
        ELSE
          IF (HAVE_DBCDY) THEN
            CALL NEWITER(NEQN,NSUB,G_X,G_Y,GE_INFO,G_K_DISCRETE,             &
                         FSUB,BCSUB,DUMMY_DF,DBCDY)
          ELSE
            CALL NEWITER(NEQN,NSUB,G_X,G_Y,GE_INFO,G_K_DISCRETE,             &
                         FSUB,BCSUB,DUMMY_DF,DUMMY_DBC)
          END IF
        END IF

        IF (GE_INFO /= 0) THEN
          IF (G_STOP_ON_FAIL) THEN
            PRINT *,'Computation of YERROR failed--the solution SOL is in doubt.'
            STOP
          END IF
        END IF
        
!       Compare discrete solution, G_Y, from doubled mesh with the original 
!       solution at points of original mesh, to estimate the error there.
!       Allocate temporary storage, ORIGY, for solution values from BVP_EVAL
!       and WT for the weights in measuring relative error.
        YERROR = 0D0
        ALLOCATE(ORIGY(NEQN),WT(NEQN),STAT=IER)
        CALL CHECK_STAT(IER)
        DO I = 1,SOL%NPTS
!         Solution values.
          CALL BVP_EVAL(SOL,SOL%X(I),ORIGY(1:NODE))
!         Parameter values.
          IF (NPAR > 0) THEN
            CALL BVP_EVAL(SOL,ORIGY(NODE+1:NEQN))
          END IF
          WT = 1D0 + ABS(ORIGY(1:NEQN))  
          YERROR = MAX(YERROR,MAXVAL(ABS(ORIGY(1:NEQN)-                    &
                          G_Y((I-1)*2*NEQN+1:(I-1)*2*NEQN+NEQN))/WT))
        END DO
        DEALLOCATE(ORIGY,STAT=IER)
        CALL CHECK_STAT(IER)

        YERROR = YERROR*(2**BVP_METHOD)/(2**BVP_METHOD - 1D0)

      END IF ! End of optional computation of YERROR
      
    END IF ! (INFO >= 0)
    
!   Deallocate global arrays before return.
    DEALLOCATE(G_X,G_Y,G_K_DISCRETE,G_K_INTERP,C_C,C_V,C_B,C_X,C_C_STAR,         &
               C_V_STAR,C_X_STAR,STAT=IER)
    CALL CHECK_STAT(IER)
    IF (SINGULAR) THEN
      DEALLOCATE(MSING,FSING,BCSING,STAT=IER)
      CALL CHECK_STAT(IER)
    END IF

  END FUNCTION BVP_SOLVER


  SUBROUTINE MIRKDC(INFO,FSUB,GSUB,DFSUB,DGSUB)
! This routine controls the overall computation. Each iteration consists of three
! main steps: the setup and solution of the nonlinear system arising from the
! boundary conditions and the discretization of the ODEs for the current mesh,
! the setup of the interpolant and the evaluation of the defect, and determination 
! of a new mesh. Success or failure is indicated by INFO.  It is 0 for a successful 
! termination, i.e., the computed solution has a defect within the user-specified 
! tolerance.  It is -1 for an unsuccessful termination; the mesh selection algorithm 
! needed to use more than the maximum number of subintervals allowed; thus the 
! solver has not been able to compute a solution that passes the test on the defect. 
! Internally, INFO indicates difficulties with the computation such as failed Newton 
! iterations, singular Newton matrices, etc. with values of 1,2,3 or 4. These are 
! recoverable error situations since the code will restart the computation on a 
! different mesh, if possible. FSUB,GSUB,DFSUB,DGSUB are user-supplied subroutines 
! for the ODEs and boundary conditions and their derivatives.
!--------------------------------------------------------------------------------
!     CALLED BY: BVP_SOLVER
!     CALLS: NEWITER,DEFECT_ESTIMATE,MESH_SELECTOR,INTERP_EVAL,HALF_MESH,
!            CHECK_STAT
!---------------------------------------------------------------------------------
!   Output arguments:
    INTEGER :: INFO

!   User-supplied subroutines:
    EXTERNAL FSUB,GSUB,DFSUB,DGSUB 
    
!   Local variables:
    INTEGER :: I ! Loop index.
    INTEGER :: IER ! Error status indicator for storage allocation.
    INTEGER :: NSUB_STAR ! Number of subintervals for new mesh.
    DOUBLE PRECISION :: DEFECT_NORM ! Weighted maximum norm of defect estimate.
    DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: Y_OLD ! Copy of current discrete 
                                                       ! solution. 
    DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: DEFECT ! Estimate of maximum 
                      ! defect for each solution component on each subinterval.
    DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: MESH_NEW ! New mesh.  
    DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: Y_NEW ! Iniital guess for  
                                                       ! solution on new mesh.                 
                                                         
!   Global variables:
!   MXNSUB is the maximum number of subintervals allowed.
!   PROFILE gives the amount of trace information about the computation that is
!     desired by the user.
!   BVP_METHOD identifies which continuous Runge-Kutta method is to be set up. 
!   BVP_TOL is input TOL or the default value.
!   NODE is the number of ODEs.
!   NPAR is the number of unknown parameters.
!   NEQN equals NODE + NPAR.
!   NSUB is the number of subintervals of the current mesh.
!   G_X, current mesh; its dimension is 0:NSUB.
!   G_Y, current discrete solution approximation, in column form; 
!     its dimension is NEQN*(NSUB+1).
!   G_K_DISCRETE contains the discrete Runge-Kutta stages.  It has dimension
!     C_S*NEQN*NSUB.
!   G_K_INTERP contains the discrete continuous Runge-Kutta stages. It has 
!     dimension (C_S_STAR-C_S)*NEQN*NSUB.
!   C_S and C_C_STAR are the number of discrete Runge-Kutta stages and 
!     continuous Runge-Kutta stages, respectively. 
!   G_STOP_ON_FAIL determines whether or not the code stops itself after a 
!     failure or simply returns with INFO=-1.
!----------------------------------------------------------------------------
    IF (PROFILE > 0) THEN
      PRINT *,' '
      PRINT *,'The order of the Runge-Kutta method is ',BVP_METHOD,'.'
      PRINT *,'The number of ODEs is ',NODE,'.'
      IF (NPAR > 0) THEN
        PRINT *,'The number of unknown parameters is ',NPAR,'.'
      END IF
      PRINT *,' '
    END IF

    IF (PROFILE > 0) THEN
      PRINT *,'The initial mesh has ',NSUB,' subintervals.' 
      PRINT *,' '
    END IF

    IF (PROFILE > 1) THEN
      WRITE(*,"(1X,A19,I5,A14)")'The initial mesh is:'
      WRITE(*,"(7F10.6)")(G_X(I),I=0,NSUB)
      PRINT *,' '
    END IF


!   INFO is already zero from initialization by BVP_INIT; 
!   initialize DEFECT_NORM so that loop will execute.
    DEFECT_NORM = 2D0*BVP_TOL
    
!   Execute loop until we have an acceptable solution or until
!   we exceed the number allowed for subintervals.
    DO WHILE ((INFO == 0) .AND. (DEFECT_NORM > BVP_TOL)) 
      IF (PROFILE > 0) THEN
          PRINT *,'Begin the Newton iteration:' 
          PRINT *,' '
      END IF

!     Allocate storage for Y_OLD, DEFECT, MESH_NEW, and Y_NEW. 
      ALLOCATE(Y_OLD(NEQN*(NSUB+1)),DEFECT(NEQN*NSUB),MESH_NEW(0:4*NSUB),        &
               Y_NEW(NEQN*(4*NSUB+1)),STAT=IER)
      CALL CHECK_STAT(IER)

!     Save a copy of the current solution approximation.
      Y_OLD = G_Y

!     Setup and solve the Newton system for the current mesh.
      CALL NEWITER(NEQN,NSUB,G_X,G_Y,INFO,G_K_DISCRETE,FSUB,GSUB,DFSUB,DGSUB)

      IF ((INFO /= 0) .AND. (PROFILE == 1)) THEN 
              WRITE(*,*)'The Newton iteration failed.'
              PRINT *,' '
      END IF
      IF ((INFO /= 0 ) .AND. (PROFILE > 1)) THEN
          IF (INFO == 1) THEN
              PRINT *,'The Newton iteration failed because the',                 &
                      ' maximum number of Newton iterations was exceeded.'
              PRINT *,' '
          END IF
          IF (INFO == 2) THEN
              PRINT *,'The Newton iteration failed because a ',                  &
                      'singular Newton matrix was encountered.'
              PRINT *,' '
         END IF
          IF (INFO == 3) THEN
              PRINT *,'The Newton iteration has been terminated because ',       &
                      'it is not converging at a satisfactory rate.' 
              PRINT *,' '
          END IF
      END IF    

!     If the Newton iteration has converged, compute the defect estimate.
      IF (INFO == 0) THEN
        CALL DEFECT_ESTIMATE(NEQN,NSUB,G_X,G_Y,DEFECT,DEFECT_NORM,INFO,          &
                             G_K_DISCRETE,G_K_INTERP,FSUB)

!       DEFECT_ESTIMATE returns the defect estimate and its max norm. If 
!       the norm is too large, set INFO to 4 to indicate that the solution 
!       should not be trusted.
        IF (PROFILE > 0) THEN
          PRINT *,'The Newton iteration was successful. '
          PRINT *,'Build a continuous approximation to the solution and',        &
                  ' sample the defect.'
          PRINT *,'The norm of the defect is ',DEFECT_NORM
          PRINT *,' '
          IF (INFO == 4) THEN
            PRINT *,'Since the defect is greater than ',                         &
                    '10%, the solution is not acceptable.'  
            PRINT *,' '
          END IF ! (INFO == 4)
        END IF ! (PROFILE > 0)          
      END IF ! (INFO == 0)

!     If the Newton iteration failed or the defect was not acceptable, INFO is 
!     nonzero and we will proceed to the ELSE clause of the following 
!     IF-THEN-ELSE statement where mesh halving will be attempted. Otherwise, 
!     we will proceed to the THEN clause where mesh redistribution/refinement 
!     will be attempted, unless the defect satisfies the tolerance.
      IF (INFO == 0) THEN                                              
!       Check user tolerance. If it is not satisfied, try again on a new mesh, 
!       if possible.
        IF (DEFECT_NORM > BVP_TOL) THEN

          IF (PROFILE > 0) THEN
            PRINT *,'User defined tolerance',BVP_TOL,' has not been satisfied.'
            PRINT *,'Construct a new mesh which equidistributes the defect.'
            PRINT *,' '
          END IF

!         Attempt to select a new mesh.
          CALL MESH_SELECTOR(NEQN,NSUB,G_X,DEFECT,NSUB_STAR,MESH_NEW,INFO)

!         If we were successful in obtaining a new mesh, compute a new estimate
!         for the solution on this mesh and then update NSUB, MESH, and Y.
          IF (INFO == 0) THEN
            IF (PROFILE > 0) THEN
              PRINT *,'The new mesh will have ',NSUB_STAR,' subintervals.'
              PRINT *,' '                              
            END IF
            IF (PROFILE > 1) THEN
              WRITE(*,"(7F10.6)")(MESH_NEW(I), I=0,NSUB_STAR)
              PRINT *,' '
            END IF

!           Use current computed solution to generate next initial guess. 
            DO I = 0, NSUB_STAR
              CALL INTERP_EVAL(NEQN,NSUB,G_X,G_Y,MESH_NEW(I),                    &
                               Y_NEW(I*NEQN+1:I*NEQN+NEQN),G_K_DISCRETE,         &
                               G_K_INTERP)
            END DO

!           Deallocate current storage for mesh and the solution and allocate
!           new storage of appropriate size.
            DEALLOCATE(G_X,G_Y,STAT=IER)
            CALL CHECK_STAT(IER)
            ALLOCATE(G_X(0:NSUB_STAR),G_Y(NEQN*(NSUB_STAR+1)),STAT=IER)
            CALL CHECK_STAT(IER)

!           Copy new mesh to G_X, new solution approximation to G_Y, and 
!           update NSUB.
            G_X(0:NSUB_STAR) = MESH_NEW(0:NSUB_STAR)
            G_Y(1:(NSUB_STAR+1)*NEQN) = Y_NEW(1:(NSUB_STAR+1)*NEQN)                        
            NSUB = NSUB_STAR
            
!           Release current storage for stages and allocate new storage 
!           of appropriate size.
            DEALLOCATE(G_K_DISCRETE,G_K_INTERP,STAT=IER)
            CALL CHECK_STAT(IER)
            ALLOCATE(G_K_DISCRETE(C_S*NEQN*NSUB),G_K_INTERP((C_S_STAR-C_S)*      &
                     NEQN*NSUB), STAT=IER)
            CALL CHECK_STAT(IER)
          END IF ! (INFO == 0)
        END IF ! (DEFECT_NORM > BVP_TOL)

        IF ( (INFO == 0) .AND. (DEFECT_NORM < BVP_TOL) .AND. (PROFILE > 0) ) THEN
          PRINT *,'The user defined tolerance',BVP_TOL,' has been satisfied.'
          PRINT *,'Successful computation.'
          PRINT *,' '
        END IF

      ELSE 
!       The Newton iteration failed or the defect was not acceptable.
        IF (PROFILE > 0) THEN
            PRINT *,'Cannot obtain a solution for the current mesh.'
            PRINT *,' '
        END IF

!       Halve the current mesh, if possible.
        IF (2*NSUB > MXNSUB) THEN
!         New mesh would be too large. Terminate the computation.
          IF (PROFILE > 0) THEN
            PRINT *,'Number of subintervals needed for new mesh would ',         &
                    'exceed current allowed maximum.'
            PRINT *,' '
          END IF
          INFO = -1
        ELSE
          IF (PROFILE > 0) THEN
            PRINT *,'Y will be reset to previous solution approximation.'
            PRINT *,' '
          END IF

!         Deallocate current storage for the solution and allocate new storage
!         of appropriate size.
          DEALLOCATE(G_Y,STAT=IER)
          CALL CHECK_STAT(IER)
          ALLOCATE(G_Y(NEQN*(2*NSUB+1)),STAT=IER)
          CALL CHECK_STAT(IER)
          
!         Reset to previous solution approximation saved in Y_OLD. The mesh has
!         been doubled, so copy the values into every other location in G_Y and
!         average the values at midpoints.
          G_Y(1:NEQN) = Y_OLD(1:NEQN)
          DO I = 1,NSUB
            G_Y(2*I*NEQN+1:2*I*NEQN+NEQN) = Y_OLD(I*NEQN+1:I*NEQN+NEQN)
            G_Y(2*I*NEQN-NEQN+1:2*I*NEQN) = (Y_OLD((I-1)*NEQN+1:(I-1)*NEQN+NEQN) &
                                             + Y_OLD(I*NEQN+1:I*NEQN+NEQN))/2D0
          END DO

!         Double number of subintervals and call HALF_MESH to compute new mesh.
          NSUB_STAR = 2*NSUB
          CALL HALF_MESH(NSUB,MESH_NEW)

!         Deallocate current storage for mesh and allocate new storage
!         of appropriate size.
          DEALLOCATE(G_X,STAT=IER)
          CALL CHECK_STAT(IER)
          ALLOCATE(G_X(0:NSUB_STAR),STAT=IER)
          CALL CHECK_STAT(IER)

!         Copy new mesh to G_X and update NSUB.
          G_X(0:NSUB_STAR) = MESH_NEW(0:NSUB_STAR) 
          NSUB = NSUB_STAR

!         Release current storage for stages and allocate new storage 
!         of appropriate size.
          DEALLOCATE(G_K_DISCRETE,G_K_INTERP,STAT=IER)
          CALL CHECK_STAT(IER)
          ALLOCATE(G_K_DISCRETE(C_S*NEQN*NSUB), &
                   G_K_INTERP((C_S_STAR-C_S)*NEQN*NSUB),STAT=IER)
          CALL CHECK_STAT(IER)

          IF (PROFILE > 0) THEN
            PRINT *,'Halve each subinterval of the current mesh and try again.'
            PRINT *,'The number of subintervals in the new mesh will be ',       & 
            NSUB,'.'
            PRINT *,' '
          END IF
          IF (PROFILE > 1) THEN
            WRITE(*,"(7F10.6)")(G_X(I), I = 0,NSUB)
            PRINT *,' '
          END IF 

!         Reset INFO and DEFECT_NORM > BVP_TOL to allow another iteration. 
          INFO = 0            
          DEFECT_NORM = 2D0*BVP_TOL

        END IF ! (2*NSUB > MXNSUB)
      END IF ! (INFO == 0)

      DEALLOCATE(Y_OLD, DEFECT, MESH_NEW, Y_NEW, STAT=IER)
      CALL CHECK_STAT(IER)

    END DO !  WHILE ((INFO == 0) .AND. (DEFECT_NORM > BVP_TOL))

    IF (INFO /= 0) THEN
!     The computation has failed. A brief report and a STOP follow unless the 
!     user has suppressed this action by setting STOP_ON_FAIL = .FALSE.
      IF (G_STOP_ON_FAIL) THEN
        PRINT *,' '
        PRINT *,'The computation has failed.'
        PRINT *,' '
        PRINT *,'1) If the profiling information begins with and simply proceeds'  
        PRINT *,'through a series of failed Newton iterations each followed by a'
        PRINT *,'doubling of the mesh, then there may be an error in the coding'
        PRINT *,'of the problem. It is recommended that the routines FSUB, GSUB,'
        PRINT *,'DFSUB, and DGSUB be checked. A better initial guess may also help.'
        PRINT *,' ' 
        PRINT *,'2) If the profiling information shows a series of successful,'
        PRINT *,'converged Newton iterations with a halt simply because the maximum'
        PRINT *,'number of subintervals was exceeded, then it is recommended that'
        PRINT *,'the program be run again with a larger value for this maximum.'
        PRINT *,'This can be done by making a call, prior to the call to BVP_SOLVER,'
        PRINT *,'to the routine BVP_INIT, with an additional argument,'
        PRINT *,'MAX_NUM_SUBINTERVALS, set equal to the desired maximum number of'
        PRINT *,'subintervals.'
        PRINT *,' '
        PRINT *,'3) More information about the computation may be obtained by choosing '
        PRINT *,'a larger value for the input argument TRACE. The maximum value for this'
        PRINT *,'input argument is 2.'

!       Since computation has failed, terminate now.
        STOP
      END IF ! (G_STOP_ON_FAIL)
      
    END IF ! (INFO \= 0)

  END SUBROUTINE MIRKDC


  SUBROUTINE CRITERION(NEQN,NSUB,MESH,Y,TOP,BLOCKS,BOT,PIVOT,PHI,                &
                       DELTA,G,K_DISCRETE,FSUB,GSUB)
! This routine evaluates the natural criterion function for the current
! iterate Y, specifically G = 0.5 * || inverse(J)*PHI(Y) ||**2.  Here
! PHI(Y) is the residual function evaluated at Y and J is the Jacobian
! of PHI evaluated at Y, or possibly at a previous value of the iterate
! if a Fixed Jacobian iteration has been employed.  Details are found in
! Ascher, Mattheij, and Russell, p. 335.  Some intermediate quantities
! are returned for use in the next step of the iteration, namely
! PHI = PHI(Y) and DELTA = inverse(J)*PHI(Y).  The input arrays 
! TOP, BLOCKS, BOT, PIVOT contain J in factored form. NEQN is the
! number of ODEs plus the number of unknown parameters; NSUB is the
! number of subintervals in the mesh, stored in MESH. By-products of
! this computation are the Runge-Kutta stages, returned from this
! routine in K_DISCRETE. FSUB and GSUB are user-supplied subroutines
! for the ODEs and boundary conditions.
!-----------------------------------------------------------------------
!     CALLED BY: FIXED_JACOB,DAMP_FACTOR
!     CALLS: RESID,CRSLVE
!-----------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NEQN,NSUB,PIVOT(NEQN*(NSUB+1))
    DOUBLE PRECISION :: MESH(0:NSUB),Y(NEQN*(NSUB+1))
    DOUBLE PRECISION :: TOP(NEQN*LEFTBC),BLOCKS(2*NEQN**2 * NSUB),               &
                        BOT(NEQN*RIGHTBC)

!   Output arguments:
    DOUBLE PRECISION :: PHI(NEQN*(NSUB+1)),DELTA(NEQN*(NSUB+1)),                 &
                        G,K_DISCRETE(MXS*NEQN*NSUB)

!   User-supplied subroutines:
    EXTERNAL FSUB,GSUB

!   Local variable:
    DOUBLE PRECISION :: PHI_COPY(NEQN*(NSUB+1))

!   Global variables:
!   MXS is the (integer) maximum number of stages in the Runge-Kutta formula.
!   LEFTBC is the number of left boundary conditions.
!   RIGHTBC is the number of right boundary conditions.
!   PROFILE specifies the level of output.
    EXTERNAL CRSLVE
!----------------------------------------------------------------------------  
!   Compute the residual PHI(Y).
    CALL RESID(NEQN,NSUB,MESH,Y,PHI,K_DISCRETE,FSUB,GSUB)

!   Compute DELTA = inverse(J)*PHI(Y) by solving a linear system.  The factored
!   matrix J is stored in TOP, BLOCKS, BOT, and PIVOT.  CRSLVE overwrites the
!   right hand side in solving the system, so we call it with a copy of PHI.
    PHI_COPY = PHI
    CALL CRSLVE(TOP,LEFTBC,NEQN,BLOCKS,NEQN,2*NEQN,NSUB,BOT,RIGHTBC,PIVOT, &
                PHI_COPY,DELTA)

    G = MAXVAL(ABS(DELTA))
    IF (G < SQRT(HUGE(1D0))) THEN
      G = 0.5D0*G**2
    ELSE
!     G is much too large.  Set a negative value for G which will be 
!     trapped in the calling routine and the iteration terminated.
      G = -1D0
      IF (PROFILE > 1) THEN
        PRINT *,'Computed Newton correction is much too large.',&
                ' Terminate current Newton interation.'
        PRINT *,' '
      END IF
    END IF

  END SUBROUTINE CRITERION


  SUBROUTINE DAMP_FACTOR(NEQN,NSUB,MESH,Y,DELTA_0,G_0,TOP,BOT,BLOCKS,PIVOT,      &
                   LAMBDA,PHI,DELTA,G,FIXED_JACOBIAN,INFO,K_DISCRETE,FSUB,GSUB)
! This routine computes a damping factor for a damped Newton step and updates
! the Newton iterate Y, which represents a discrete solution to the BVP on a 
! mesh, MESH, of NSUB subintervals. NEQN is the number of ODEs plus the number 
! of unknown parameters. On input LAMBDA is an estimate for the damping factor; 
! on output, it contains the damping factor used to perform the update. TOP, BOT, 
! BLOCKS, and PIVOT provide the Newton matrix in factored form. DELTA_0 is the 
! Newton correction for the input Y. G_0 is the natural criterion function value
! corresponding to the input Y. PHI is the residual function value for the new 
! iterate. DELTA is the Newton correction for the new iterate and G is the 
! corresponding natural criterion function value. On output FIXED_JACOBIAN is 
! .TRUE. if the next step should be a fixed Jacobian step and .FALSE. if the next 
! step should be another damped Newton step. Success is indicated by INFO = 0. 
! A large natural criterion function value is indicated by INFO = 3. If it is 
! impossible to find a suitable damping factor, we set INFO = 3. A by-product of 
! this computation is the Runge-Kutta stages, which are returned in K_DISCRETE. 
! FSUB and GSUB are provided by the user to define the ODEs and boundary 
! conditions.
!------------------------------------------------------------------------
!     CALLED BY: DAMPED_NEWT
!     CALLS: CRITERION
!------------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NEQN,NSUB, PIVOT(NEQN*(NSUB+1))  
    DOUBLE PRECISION :: MESH(0:NSUB),Y(NEQN*(NSUB+1)),DELTA_0(NEQN*(NSUB+1)),    &
                         G_0,TOP(NEQN**2), BOT(NEQN**2),BLOCKS(2*NEQN**2*NSUB),  &
                         LAMBDA
!   Output arguments:
    LOGICAL :: FIXED_JACOBIAN
    INTEGER :: INFO
    DOUBLE PRECISION :: PHI(NEQN*(NSUB+1)),DELTA(NEQN*(NSUB+1)),G,               &
                        K_DISCRETE(MXS*NEQN*NSUB)

!   User-supplied subroutines:
    EXTERNAL FSUB,GSUB
 
!   Local variables:
    LOGICAL :: ACCEPT ! Used to monitor search for new damping factor.
    DOUBLE PRECISION, PARAMETER :: SIGMA=0.01D0,TAU=0.1D0 
!   SIGMA ensures that reduction in size of natural criterion function value 
!   for new iterate will be nonnegligible. TAU controls how much values for 
!   LAMBDA are allowed to change from one step to the next in the iteration 
!   performed by this routine. 
    DOUBLE PRECISION :: Y_0(NEQN*(NSUB+1)) ! Current iterate saved here while 
                                           ! new iterate is determined.
!   Global variables:
!   MXS is the maximum number of Runge-Kutta stages.
!   PROFILE gives the level of output to be provided.
!   LAMBDA_MIN is the minimum value for the damping factor.
!   SINGULAR specifies whether or not the BVP has a singular term.
!   BCSING is a quantity associated with the treatment of a singular term.
!   NODE is the number of ODEs.
!--------------------------------------------------------------------------
!   Iterative determination of new damping factor.
    Y_0 = Y
    ACCEPT = .FALSE.
    DO WHILE ( (.NOT. ACCEPT) .AND. (INFO == 0) )
!     Update iterate value using input damping factor and Newton correction.
      Y = Y_0 - LAMBDA*DELTA_0

      IF (SINGULAR) THEN
        ! Impose necessary condition on approximation at left end point.
        Y(1:NODE) = MATMUL(BCSING,Y(1:NODE))
      END IF
             
!     Measure suitability of current damping factor by seeing how much the
!     trial iterate reduces the natural criterion function. CRITERION returns 
!     G, PHI, and DELTA = inv(J)*PHI, with J evaluated at the previous iterate.
      CALL CRITERION(NEQN,NSUB,MESH,Y,TOP,BLOCKS,BOT,PIVOT,PHI,DELTA,G,  &
                     K_DISCRETE,FSUB,GSUB)

!     If natural criterion function value has overflowed, CRITERION will 
!     trap this and return G = -1. Set INFO = 3 here to cause a return.
      IF (G < 0D0) THEN
        INFO = 3
      ELSE
!       Check G to ensure that reduction in size is sufficient. If so, LAMBDA
!       and thus Y are acceptable. 
        IF (G <= (1D0 - 2D0*LAMBDA*SIGMA)*G_0) THEN
          ACCEPT = .TRUE.
          IF ( ABS(LAMBDA - 1D0) <= 0D0 ) FIXED_JACOBIAN = .TRUE.
        ELSE
!         Reduction in G is not sufficient. Compute a new damping factor.
          LAMBDA = MAX(TAU*LAMBDA,(LAMBDA**2 * G_0)/((2*LAMBDA-1D0)*G_0+G))
          IF (PROFILE > 1) THEN
            PRINT *,'Damped Newton step: value of damping factor is ', LAMBDA
            PRINT *,' '
          END IF

!         If new LAMBDA is too small, set INFO = 3 to cause a return.
          IF (LAMBDA < LAMBDA_MIN) INFO = 3
        END IF ! (G <= (1D0-2D0*LAMBDA*SIGMA)*G_0)
      END IF ! (G < 0D0)
    END DO

  END SUBROUTINE DAMP_FACTOR

  SUBROUTINE DAMPED_NEWT(NEQN,NSUB,MESH,Y,LAMBDA,PHI,TOP,BOT,BLOCKS,PIVOT,&
                         FIXED_JACOBIAN,CONVRG,DELTA,DELTA_0_NORM,G,INFO,&
                         K_DISCRETE,FSUB,GSUB,DFSUB,DGSUB)
! This routine takes one damped Newton step. Y contains the current iterate 
! and PHI contains the corresponding residual function value for the current 
! mesh of NSUB subintervals stored in MESH. NEQN is the number of ODEs plus 
! the number of unknown parameters. LAMBDA is the damping factor. TOP, BOT, 
! BLOCKS, and PIVOT, store the Newton matrix in factored form. FIXED_JACOBIAN 
! indicates whether the next step should be a fixed Jacobian step and CONVRG 
! indicates whether the iteration has converged. DELTA is the Newton correction 
! vector and DELTA_0_NORM is its norm. G is the value of natural criterion 
! function. A successful return is indicated by INFO = 0. If a singular matrix 
! is encountered, INFO = 2. A large natural criterion function value is indicated
! by INFO = 3. A by-product of this computation is the Runge-Kutta stages which 
! are returned in the array K_DISCRETE. The subroutines FSUB, GSUB, DFSUB, and 
! DGSUB are provided by the user to define the ODEs and boundary conditions and 
! their partial derivatives.
!------------------------------------------------------------------------
!     CALLED BY: NEWITER
!     CALLS: NEWMAT,COLROW,DAMP_FACTOR
!----------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NEQN,NSUB
    DOUBLE PRECISION :: MESH(0:NSUB),Y(NEQN*(NSUB+1)),LAMBDA,                    &
                        PHI(NEQN*(NSUB+1)),DELTA(NEQN*(NSUB+1)),DELTA_0_NORM

!   Output arguments:
    INTEGER :: INFO,PIVOT(NEQN*(NSUB+1))
    DOUBLE PRECISION :: TOP(NEQN**2),BOT(NEQN**2),BLOCKS(2*NEQN**2*NSUB),        &
                        G,K_DISCRETE(MXS*NEQN*NSUB)
    LOGICAL :: FIXED_JACOBIAN, CONVRG
      
!   User-supplied subroutines:
    EXTERNAL FSUB,GSUB,DFSUB,DGSUB

!   Local variables:
    INTEGER :: FACTOR ! Trouble flag for call to COLROW.
    DOUBLE PRECISION :: G_0 ! Natural criterion function value for input Y.
    DOUBLE PRECISION :: MU ! Predicted value for new damping factor.
    DOUBLE PRECISION :: DELTA_0(NEQN*(NSUB+1)) ! Copy of Newton correction 
                                               ! vector.
      
!   Global variables:
!   MXS is the number of maximum number of Runge-Kutta stages.
!   LEFTBC is the number of left boundary conditions.
!   RIGHTBC is the number of right boundary conditions.
!   PROFILE defines the level of output to be provided.
!   NEWTON_TOL is the tolerance on the Newton iteration.
!   LAMBDA_MIN is the minimum value for the damping factor.
    EXTERNAL COLROW
!----------------------------------------------------------------------------
    INFO = 0  

!   Form the Newton matrix with NEWMAT.
    CALL NEWMAT(LEFTBC,RIGHTBC,NEQN,NSUB,MESH,Y,TOP,BLOCKS,BOT, &
                K_DISCRETE,FSUB,GSUB,DFSUB,DGSUB)
                
!   Factor the Newton matrix and solve a linear system with the input
!   vector PHI as right hand side to get the Newton correction DELTA_0.
    CALL COLROW((NSUB+1)*NEQN,TOP,LEFTBC,NEQN,BLOCKS, &
                NEQN,2*NEQN,NSUB,BOT,RIGHTBC,PIVOT,PHI,DELTA_0,FACTOR)
    IF ( FACTOR == -1 ) THEN
      INFO = 2 ! Singular Newton matrix encountered.
    END IF

    IF (INFO == 0) THEN
!     Compute natural criterion function value.
      G_0 = MAXVAL(ABS(DELTA_0))
      IF (G_0 > SQRT(HUGE(1D0))) THEN
        IF (PROFILE > 1) THEN
          PRINT *,'Computed Newton correction is much too large.',               &
                  ' Terminate current Newton interation.'
          PRINT *,' '
          INFO = 3
        END IF
      ELSE
        G_0 = 0.5D0 * G_0**2
      END IF
    END IF
      
    IF (INFO == 0) THEN
!     Test if norm of scaled Newton correction, DELTA_0, is sufficiently small
!     to signal convergence. 
      CONVRG = ALL( ABS(DELTA_0) <= NEWTON_TOL*(ABS(Y) + 1D0) )
      IF (CONVRG) THEN
!       Update iterate with full Newton correction.
        Y = Y - DELTA_0
      ELSE           
!       Attempt to update the Newton iterate with a damped Newton correction. 
!       Compute a prediction for new LAMBDA based on previous LAMBDA value. 
!       If the input LAMBDA = 0, the previous iteration was not a damped Newton 
!       step and appropriate information for the prediction of a new LAMBDA 
!       value is not available.
        IF ( ABS(LAMBDA) <= 0D0 ) THEN
          LAMBDA = 1D0
        ELSE
!         Predict LAMBDA based on previous LAMBDA, DELTA, and DELTA_0_NORM.
          DELTA = DELTA - DELTA_0
          MU = (LAMBDA*DELTA_0_NORM) / MAXVAL(ABS(DELTA))
          LAMBDA = MAX(LAMBDA_MIN,MIN(1D0,MU))
        END IF
        IF (PROFILE > 1) THEN
          PRINT *,'Damped Newton iteration: value of initial damping',           &
                  ' factor is', LAMBDA
          PRINT *,' '
        END IF
        
!       Update norm of Newton correction.
        DELTA_0_NORM = SQRT(2D0*G_0)

!       Compute a new iterate by updating Y with a damped Newton correction.
        CALL DAMP_FACTOR(NEQN,NSUB,MESH,Y,DELTA_0,G_0,TOP,BOT,BLOCKS,PIVOT,      & 
                         LAMBDA,PHI,DELTA,G,FIXED_JACOBIAN,INFO,K_DISCRETE,      &
                         FSUB,GSUB)
                         
!       Upon return DAMP_FACTOR has computed an appropriate damping factor, 
!       LAMBDA, and updated the iterate, Y, accordingly. This routine has also
!       updated the corresponding residual function value, PHI, the Newton 
!       correction DELTA, and the natural criterion function value, G, all for 
!       the newly updated Y.

        IF ((INFO == 0) .AND. (PROFILE > 1)) THEN 
          PRINT *,'Norm of damped Newton correction is',LAMBDA*DELTA_0_NORM
          PRINT *,' '
        END IF
      END IF ! IF (CONVERG)
    END IF ! IF (INFO==0)

  END SUBROUTINE DAMPED_NEWT


  SUBROUTINE DEFECT_ESTIMATE(NEQN,NSUB,MESH,Y,DEFECT,DEFECT_NORM,                &
                             INFO,K_DISCRETE,K_INTERP,FSUB) 
! This routine uses the discrete approximate solution Y on the mesh stored in
! MESH of NSUB subintervals together with the Runge-Kutta stages stored
! in K_DISCRETE and some new stages stored in K_INTERP to construct a
! continuous approximate solution on the whole interval. NEQN is the number 
! of differential equations plus the number of unknown parameters. The defect
! of the continuous extension is estimated on each subinterval for each solution 
! component. The Ith block of NEQN locations in the vector DEFECT holds the 
! defect information associated with subinterval I. A weighted maximum norm of 
! DEFECT is returned in DEFECT_NORM. A successful return is indicated by  
! INFO = 0. If DEFECT_NORM is bigger than 10%, INFO is set to 4 to cause  
! rejection of the current Y. The subroutine FSUB is provided by the user 
! to define the ODEs.
!
! The defect on subinterval I evaluated at TAU is 
!    Z'(MESH(I-1)+TAU*HI) - F(MESH(I-1)+TAU*HI,Z(MESH(I-1)+TAU*HI)).
! A weighted measure of the defect at this point, namely |DEFECT|/(|F| + 1), 
! is computed. The maximum defect on each subinterval is estimated by the 
! larger of the defect values at TAU as just described and at 1-TAU.  The value 
! of TAU is set as a global variable appropriate to the method.
!
! The boundary conditions are solved during the computation of the discrete 
! solution. Thus the boundary conditions will be satisfied to within the Newton 
! tolerance, which is well below the user's tolerance.  Accordingly, the defect 
! of the boundary conditions will be less than the user's tolerance for all 
! converged discrete solutions, and there is no need to monitor it here. 
!-----------------------------------------------------------------------------
!     CALLED BY: MIRKDC
!     CALLS: INTERP_WEIGHTS,INTERP_SETUP,SUM_STAGES,P_FSUB
!------------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NEQN,NSUB
    DOUBLE PRECISION :: MESH(0:NSUB),Y(NEQN*(NSUB+1)),K_DISCRETE(MXS*NEQN*NSUB)

!   Output arguments:
    INTEGER :: INFO
    DOUBLE PRECISION :: DEFECT(NSUB*NEQN),DEFECT_NORM,K_INTERP(MXS*NEQN*NSUB)
    
!   User-supplied subroutines:
    EXTERNAL FSUB

!   Local variables:
    INTEGER :: I,J ! Loop indices.
    DOUBLE PRECISION :: HI ! Size of the i-th subinterval.
    DOUBLE PRECISION :: WEIGHTS(MXS,MXS) ! Weight polynomials evaluated at
                                         ! sample points.
    DOUBLE PRECISION :: WEIGHTS_PRIME(MXS,MXS) ! Derivatives of weight 
                                               ! polynomials evaluated at 
                                               ! sample points.
    DOUBLE PRECISION MAX_ESTIMATE, ESTIMATE(MXS) ! Maximum defect and estimates.                
    DOUBLE PRECISION :: Z(NEQN) ! Argument for f evaluation.
    DOUBLE PRECISION :: Z_PRIME(NEQN) ! Derivative of interpolant.
    DOUBLE PRECISION :: F_I(NEQN),DEFECT_I(NEQN) ! Temporary storage.
    DOUBLE PRECISION :: SAMPLE(2) ! Sample points.

!   Global variables:
!   MXS is the maximum number of stages of the Runge-Kutta method.
!   C_S is the number of discrete Runge-Kutta stages. 
!   C_S_STAR is the total number of stages required. It is the number of discrete
!   stages plus the number of additional stages required for the interpolant. 
!   C_TAU_STAR  is the relative position of one of the sample points within each 
!   subinterval; the other sample point has the relative position 1-C_TAU_STAR.
!--------------------------------------------------------------------------------
!   Setup information about sample points.
    SAMPLE = (/ C_TAU_STAR, 1D0-C_TAU_STAR /)
      
!   Setup the weights for the interpolant and the derivative. 
    CALL INTERP_WEIGHTS(C_S_STAR,SAMPLE(1),WEIGHTS(:,1),WEIGHTS_PRIME(:,1))
    CALL INTERP_WEIGHTS(C_S_STAR,SAMPLE(2),WEIGHTS(:,2),WEIGHTS_PRIME(:,2))

    DO I = 1, NSUB
!     Setup new stages for this subinterval.
      HI = MESH(I) - MESH(I-1)
      CALL INTERP_SETUP(NEQN,MESH(I-1),HI,Y((I-1)*NEQN+1),Y(I*NEQN+1),C_S,       &
                        K_DISCRETE((I-1)*C_S*NEQN+1),C_S_STAR,                   &   
                        K_INTERP((I-1)*(C_S_STAR-C_S)*NEQN+1),FSUB)
                 
      MAX_ESTIMATE = 0D0
      DO J = 1,2
!       Compute value of interpolant and its first derivative at sample 
!       point J by calling SUM_STAGES. The interpolant is returned in Z;
!       the first derivative in Z_PRIME.
        CALL SUM_STAGES(NEQN,HI,Y((I-1)*NEQN+1),C_S,K_DISCRETE((I-1)*C_S*NEQN+1),&
                        C_S_STAR,K_INTERP((I-1)*(C_S_STAR-C_S)*NEQN+1),          &
                        WEIGHTS(:,J),Z,WEIGHTS_PRIME(:,J),Z_PRIME)

!       Evaluate f(Z) and return in F_I.
        CALL P_FSUB(FSUB,NEQN,MESH(I-1)+SAMPLE(J)*HI,Z,F_I)

!       Compute scaled defect and store in DEFECT_I. 
        DEFECT_I = (Z_PRIME - F_I)/(ABS(F_I) + 1D0)

!       Compute norm of defect and update the maximum defect.
        ESTIMATE(J) = MAXVAL(ABS(DEFECT_I))
        IF (ESTIMATE(J) > MAX_ESTIMATE) THEN
          MAX_ESTIMATE = ESTIMATE(J)
          DEFECT((I-1)*NEQN+1:I*NEQN) = DEFECT_I
        END IF
      END DO ! J = 1,2
    END DO ! I = 1,NSUB

    DEFECT_NORM = MAXVAL(ABS(DEFECT))

!   If the defect is too big, the current solution cannot be trusted,
!   so start over as if the Newton iteration had failed.
    IF (DEFECT_NORM > 0.1D0 ) INFO = 4
    
  END SUBROUTINE DEFECT_ESTIMATE


  SUBROUTINE FIXED_JACOB(NEQN,NSUB,MESH,Y,DELTA,G,PHI,TOP,BOT,BLOCKS,            &
                         PIVOT,FIXED_JACOBIAN,LAMBDA,CONVRG,INFO,                &
                         K_DISCRETE,FSUB,GSUB)
! This routine performs one fixed Jacobian iteration, i.e., one Newton step
! using a Newton iteration matrix evaluated at a previous iterate. Y contains 
! the current iterate and PHI contains the corresponding residual function 
! value for the current mesh, MESH, of NSUB subintervals. NEQN is the number 
! of ODEs plus the number of unknown parameters. TOP, BOT, BLOCKS, and PIVOT
! store the Newton matrix from a previous iteration in factored form. 
! FIXED_JACOBIAN indicates whether the next step should be a fixed Jacobian 
! step. CONVRG indicates whether the iteration converged. DELTA is the Newton 
! correction vector. INFO = 0 indicates a successful computation. INFO = 3
! indicates a failure as signalled by a large natural criterion function value 
! G. A by-product of this computation is the Runge-Kutta stages, which are 
! returned in the array K_DISCRETE. FSUB and GSUB define the ODEs and boundary 
! conditions.
!
! When a fixed Jacobian step is to be taken, the following quantities are
! assumed to be available from the previous iteration step: i) the current 
! iterate Y (ii) the residual function PHI evaluated at Y (iii) the Newton 
! correction vector DELTA = inverse of J *PHI, where the Jacobian J was 
! evaluated at some previous iterate (iv) the value of the natural criterion 
! function G. These quantities are provided by DAMPED_NEWT if the previous step 
! was a damped Newton step, or FIXED_JACOB if the previous step was a fixed
! Jacobian step. 
!
! The fixed Jacobian step first computes a new iterate in the local array
! Y_HAT = Y - DELTA. The acceptability of Y_HAT is determined by evaluating 
! the natural criterion function at Y_HAT.  The value G_HAT is compared to the 
! previous value of the natural criterion function. If it is sufficiently  
! smaller, the iterate is accepted and another fixed Jacobian step is attempted.  
! If not, the code returns with FIXED_JACOBIAN set to FALSE and LAMBDA set to 
! 0 to indicate that information for the prediction of a new damping factor for  
! the upcoming damped Newton step is not available.
!------------------------------------------------------------------------------
!     CALLED BY : NEWITER
!     CALLS : CRITERION,RESID
!------------------------------------------------------------------------------     
!   Input arguments:
    INTEGER :: NEQN,NSUB,PIVOT(NEQN*(NSUB+1))
    DOUBLE PRECISION :: MESH(0:NSUB),Y(NEQN*(NSUB+1)),DELTA(NEQN*(NSUB+1)),      &
                        G,PHI(NEQN*(NSUB+1)),TOP(NEQN**2),BOT(NEQN**2),          &
                        BLOCKS(2*NEQN**2*NSUB)

!   Output arguments:
    LOGICAL :: FIXED_JACOBIAN,CONVRG
    INTEGER :: INFO
    DOUBLE PRECISION :: LAMBDA,K_DISCRETE(MXS*NEQN*NSUB)

!   User-supplied subroutines:
    EXTERNAL FSUB,GSUB
! 
!   Local variables:
    DOUBLE PRECISION :: G_HAT ! value of the natural criterion function 
                              ! at the trial iterate.
    DOUBLE PRECISION :: RHO ! used in measuring performance of fixed 
                            ! jacobian iteration.  
    DOUBLE PRECISION :: Y_HAT(NEQN*(NSUB+1)) ! trial iterate value.
    DOUBLE PRECISION :: PHI_HAT(NEQN*(NSUB+1)) ! Residual function value 

!   Global variables:
!   MXS is the number of maximum number of Runge-Kutta stages.
!   PROFILE gives the level of output to be provided.
!   NEWTON_TOL is the Newton tolerance.
!------------------------------------------------------------------------------     
    INFO = 0
    RHO = 0.5D0

    IF (PROFILE > 1) THEN
      PRINT *,'Norm of Newton correction', MAXVAL(ABS(DELTA))
      PRINT *,' '
    END IF
      
!   If the Newton correction, DELTA, is less than the tolerance, the 
!   iteration has converged. 
    CONVRG = ALL( ABS(DELTA) <= NEWTON_TOL*(ABS(Y) + 1D0))

    IF (CONVRG) THEN
      Y = Y - DELTA
    ELSE
!     Compute trial iterate. 
      Y_HAT = Y - DELTA
                  
!     Compute natural criterion function to assess the trial iterate. 
!     By-products of this calculation are PHI_HAT and DELTA evaluated at Y_HAT. 
!     These MAY be useful for the next iteration step, but if Y_HAT is not 
!     acceptable, values of PHI and G associated with the iterate Y will be 
!     needed. Here we take care that they are not overwritten by providing
!     only copies Y_HAT, PHI_HAT, and G_HAT to CRITERION.
      CALL CRITERION(NEQN,NSUB,MESH,Y_HAT,TOP,BLOCKS,BOT,PIVOT,PHI_HAT,          &
                     DELTA,G_HAT,K_DISCRETE,FSUB,GSUB)

!     If the natural criterion function value is about to overflow, the CRITERION
!     subroutine will set G_HAT = -1 ( a value which cannot otherwise arise).
!     Set INFO = 3 to report failure.
      IF (G_HAT < 0D0) THEN
        INFO = 3
      ELSE
        IF ( G_HAT <= RHO * G ) THEN
!         The new iterate is acceptable.
          Y = Y_HAT
          PHI = PHI_HAT
          G = G_HAT
        ELSE
!         The new iterate is not acceptable. Signal a damped Newton step. 
          FIXED_JACOBIAN = .FALSE.
          LAMBDA = 0D0

!         If G_HAT > RHO * G, but G_HAT < G, use the new iterate and 
!         corresponding residual to begin the subsequent full Newton step. 
          IF (G_HAT < G) THEN
            Y = Y_HAT
            PHI = PHI_HAT
          ELSE
!           Begin the full Newton step using the values in Y and PHI from 
!           the beginning of this step. A call to RESID is made here to 
!           obtain the K_DISCRETE values, evaluated at Y, as these must be 
!           up-to-date with respect to Y before the damped Newton step begins.
            CALL RESID(NEQN,NSUB,MESH,Y,PHI,K_DISCRETE,FSUB,GSUB)
          END IF ! (G_HAT < G)
          
        END IF ! (G_HAT <= RHO*G)
        
      END IF ! (G_HAT < 0D0)
      
    END IF ! (CONVRG)

  END SUBROUTINE FIXED_JACOB


  SUBROUTINE HALF_MESH(NSUB,MESH_NEW)
! This routine takes as input a mesh of NSUB subintervals stored in the
! globally defined vector G_X.  It returns in MESH_NEW a new mesh which 
! has 2*NSUB subintervals obtained by halving each subinterval of G_X.
!------------------------------------------------------------------------------
!     CALLED BY: MIRKDC, MESH_SELECTOR
!------------------------------------------------------------------------------
!   Input variable:
    INTEGER :: NSUB

!   Output variable:
    DOUBLE PRECISION :: MESH_NEW(0:2*NSUB)

!   Global variable:
!   The global vector G_X is a mesh of NSUB subintervals.
!----------------------------------------------------------------------------  
    MESH_NEW(0:2*NSUB:2)   = G_X(0:NSUB)
    MESH_NEW(1:2*NSUB-1:2) = (G_X(1:NSUB) + G_X(0:NSUB-1) )/2D0

  END SUBROUTINE HALF_MESH


  SUBROUTINE INTERP_EVAL(NEQN,NSUB,MESH,Y,T,Z,K_DISCRETE,K_INTERP)
! This routine uses the current discrete solution approximation Y 
! on a mesh of NSUB subintervals stored in MESH plus additional stage 
! information associated with the underlying continuous Runge-Kutta 
! method provided in K_DISCRETE and K_INTERP to evaluate the continuous
! solution approximation at T.  The result is returned in Z, a vector of 
! length NEQN, where NEQN is the number of differential equations, NODE, 
! plus the number of unknown parameters, NPAR. The first NODE locations 
! of Z contain the approximation to the solution and the remaining locations 
! contain the approximation to the unknown parameters.
!---------------------------------------------------------------------------
!     CALLED BY: MIRKDC
!     CALLS: INTERVAL,INTERP_WEIGHTS,SUM_STAGES
!---------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NEQN,NSUB
    DOUBLE PRECISION :: MESH(0:NSUB),Y(NEQN*(NSUB+1)),T,                         &
                        K_DISCRETE(MXS*NEQN*NSUB),K_INTERP(MXS*NEQN*NSUB)
!   Output arguments:
    DOUBLE PRECISION :: Z(NEQN)

!   Local variables:
    INTEGER :: I ! Index of subinterval.
    DOUBLE PRECISION :: HI ! Size of subinterval I.
    DOUBLE PRECISION :: TAU ! Relative postion of T in subinterval I.
    DOUBLE PRECISION :: WEIGHTS(MXS) ! Weight polynomials evaluated at TAU.
  
!   Global variables:
!   MXS is the maximum number of Runge-Kutta stages.
!   C_S and C_C_STAR, are number of discrete Runge-Kutta stages and number of 
!   continuous Runge-Kutta stages, respectively. 
!-----------------------------------------------------------------------------                                           
!   The INTERVAL subroutine returns the index of an interval containing T.
    CALL INTERVAL(NSUB,MESH,T,I)
                                                                        
!   Compute TAU, the relative location of T in subinterval I.
    HI = MESH(I) - MESH(I-1)
    TAU = (T - MESH(I-1))/HI

!   Evaluate WEIGHTS at TAU.
    CALL INTERP_WEIGHTS(C_S_STAR,TAU,WEIGHTS)

!   Evaluate interpolant at TAU by taking weighted sum of Runge-Kutta stages 
!   through a call to SUM_STAGES. Result is returned in Z. 
    CALL SUM_STAGES(NEQN,HI,Y((I-1)*NEQN+1),C_S,K_DISCRETE((I-1)*C_S*NEQN+1),    &
                    C_S_STAR,K_INTERP((I-1)*(C_S_STAR-C_S)*NEQN+1),WEIGHTS,Z)

  END SUBROUTINE INTERP_EVAL

      
  SUBROUTINE INTERP_SETUP(NEQN,TIM1,HI,YIM1,YI,S,KI_DISCRETE,S_STAR,KI_INTERP,   &
                          FSUB)
! This routine constructs the S_STAR-S extra stages needed by the interpolant for 
! a single subinterval. S is the number of stages of the discrete Runge-Kutta  
! method; S_STAR is the number of stages of the continuous Runge-Kutta method. 
! NEQN is the number of differential equations plus the number of unknown  
! parameters. HI is the width of the subinterval; TIM1 is t_{i-1}, the left  
! endpoint of the subinterval; YIM1 is y_{i-1}, the solution approximation at  
! t_{i-1}; YI is y_i, the solution approximation at t_i. KI_DISCRETE contains  
! the S stages of the discrete formula Runge-Kutta on subinterval i. The 
! interpolant will use these together with the extra stages computed by this 
! routine, which are returned in KI_INTERP. FSUB is provided by the user to 
! define the ODEs.
!-------------------------------------------------------------------------------
!     CALLED BY: DEFECT_ESTIMATE
!     CALLS: P_FSUB
!---------------------------------------------------------------------------- 
!   Input arguments:
    INTEGER :: NEQN,S,S_STAR
    DOUBLE PRECISION :: TIM1,HI,YIM1(NEQN),YI(NEQN),KI_DISCRETE(S*NEQN)

!   Output arguments:
    DOUBLE PRECISION :: KI_INTERP((S_STAR-S)*NEQN)

!   User-supplied subroutines:
    EXTERNAL FSUB

!   Local variables:
    INTEGER :: R,J ! Loop indexes.
    DOUBLE PRECISION :: Y_ARG(NEQN) ! Argument for the function evaluation for 
                                    ! each new stage.
!   Global variables:
!   C_C_STAR, C_V_STAR, C_X_STAR are Runge-Kutta coefficients.
!----------------------------------------------------------------------------
    DO R = 1,S_STAR-S
      
      Y_ARG = 0D0

!     Contributions from the discrete stages.                    
      DO J = 1,S
        Y_ARG = Y_ARG + C_X_STAR(R,J)*KI_DISCRETE((J-1)*NEQN+1:J*NEQN)
      END DO

!     Contributions from the previous computed new stages.
      DO J = 1,R-1
        Y_ARG = Y_ARG + C_X_STAR(R,S+J)*KI_INTERP((J-1)*NEQN+1:J*NEQN)
      END DO

!     Complete computation of Y_ARG.
      Y_ARG = C_V_STAR(R)*YI + (1D0-C_V_STAR(R))*YIM1 + HI*Y_ARG

!     Evaluate P_FSUB at Y_ARG and TIM1+C_C_STAR(R)*HI to get new stage. 
!     Stage value is stored in KI_INTERP.
      CALL P_FSUB(FSUB,NEQN,TIM1+C_C_STAR(R)*HI,Y_ARG,KI_INTERP((R-1)*NEQN+1))
  
    END DO
    
  END SUBROUTINE INTERP_SETUP


  SUBROUTINE INTERP_TABLEAU
! This routine defines the extra coefficients for the continuous extension of the
! discrete Runge-Kutta scheme, i.e., the interpolant. The global variable 
! BVP_METHOD indicates which Runge-Kutta formula is to be used. This routine also 
! defines the sample point for the defect associated with this interpolant and 
! the order of the interpolant.
!--------------------------------------------------------------------------------
!     CALLED BY: BVP_SOLVER
!     CALLS: CHECK_STAT
!--------------------------------------------------------------------------------
!   Local variable: 
    INTEGER IER ! error flag for allocation call.
!
!   Global variables:
!   BVP_METHOD identifies which continuous Runge-Kutta method is to be set up. 
!   C_S, C_S_STAR are the number of discrete Runge-Kutta stages and the number of 
!     continuous Runge-Kutta stages, respectively. 
!   C_C_STAR, C_V_STAR, C_X_STAR are Runge-Kutta coefficients for the 
!     continuous scheme.
!   C_P is the order of the continuous Runge-Kutta formula.  
!   C_TAU_STAR is the relative sample point location for the defect 
!     sampling process. 
!--------------------------------------------------------------------------------
    ALLOCATE(C_C_STAR(C_S_STAR-C_S),C_V_STAR(C_S_STAR-C_S),                      &
             C_X_STAR(C_S_STAR-C_S,C_S_STAR),STAT=IER)
    CALL CHECK_STAT(IER)
   
    SELECT CASE (BVP_METHOD)
      CASE (2)
!       The interpolant for the 2-stage, 2nd order MIRK, stage order 1, requires  
!       no extra stages. C_S_STAR has been set to 2 in BVP_SOLVER.
        C_TAU_STAR = 0.25D0
        C_P = 2
      CASE (4)
!       Define coefficients for a 4th order, 4-stage, stage order 3 continuous 
!       scheme. There is one extra stage. C_S_STAR has been set to 4 in 
!       BVP_SOLVER.
        
        C_C_STAR(1) = 2D0/5D0
        
        C_V_STAR(1) = 2D0/5D0
        
!       Column 1 of C_X_STAR
        C_X_STAR(1,1) = 17D0/125D0

!       Column 2 of C_X_STAR
        C_X_STAR(1,2) = -13D0/125D0

!       Column 3 of C_X_STAR
        C_X_STAR(1,3) = -4D0/125D0

!       Column 4 of C_X_STAR
        C_X_STAR(1,4) = 0D0

        C_TAU_STAR = 0.226D0
        C_P = 4
      CASE (6)
!       Define coefficients for an 8-stage, 6th order, stage order 3 continuous 
!       scheme. There are three extra stages. C_S_STAR has been set to 8 in 
!       BVP_SOLVER.

        C_C_STAR(1) = 1D0/2D0
        C_C_STAR(2) = 1D0/2D0 - SQRT(7D0)/14D0
        C_C_STAR(3) = 87D0/100D0

        C_V_STAR(1) = 1D0/2D0
        C_V_STAR(2) = 1D0/2D0 - SQRT(7D0)/14D0
        C_V_STAR(3) = 87D0/100D0

!       Column 1 of C_X_STAR
        C_X_STAR(1,1) = 1D0/64D0
        C_X_STAR(2,1) = 3D0/112D0 + 9D0/1960D0*SQRT(7D0)
        C_X_STAR(3,1) = 2707592511D0/1000000000000D0 - 1006699707D0/             &
        1000000000000D0*SQRT(7D0)

!       Column 2 of C_X_STAR
        C_X_STAR(1,2) = -1D0/64D0
        C_X_STAR(2,2) = -3D0/112D0 + 9D0/1960D0*SQRT(7D0)
        C_X_STAR(3,2) = -51527976591D0/1000000000000D0 - 1006699707D0/           &
        1000000000000D0*SQRT(7D0)

!       Column 3 of C_X_STAR
        C_X_STAR(1,3) = 7D0/192D0*SQRT(21D0)
        C_X_STAR(2,3) = 11D0/840D0*SQRT(7D0) + 3D0/112D0*SQRT(7D0)*SQRT(3D0)
        C_X_STAR(3,3) = -610366393D0/75000000000D0 + 7046897949D0/               &
        1000000000000D0*SQRT(7D0) + 14508670449D0/1000000000000D0                &
        *SQRT(7D0)*SQRT(3D0)

!       Column 4 of C_X_STAR
        C_X_STAR(1,4) = -7D0/192D0*SQRT(21D0)
        C_X_STAR(2,4) = 11D0/840D0*SQRT(7D0 )- 3D0/112D0*SQRT(7D0)*SQRT(3D0)
        C_X_STAR(3,4) = -610366393D0/75000000000D0 + 7046897949D0/               &
        1000000000000D0*SQRT(7D0) - 14508670449D0/1000000000000D0                &
        *SQRT(7D0)*SQRT(3D0)

!       Column 5 of C_X_STAR
        C_X_STAR(1,5) = 0D0
        C_X_STAR(2,5) = 88D0/5145D0*SQRT(7D0)
        C_X_STAR(3,5) = -12456457D0/1171875000D0 + 1006699707D0/                 &
        109375000000D0*SQRT(7D0)

!       Column 6 of C_X_STAR
        C_X_STAR(1,6) = 0D0
        C_X_STAR(2,6) = -18D0/343D0*SQRT(7D0)
        C_X_STAR(3,6) = 3020099121D0/437500000000D0*SQRT(7D0) +                  &
        47328957D0/625000000D0

!       Column 7 of C_X_STAR
        C_X_STAR(1,7) = 0D0
        C_X_STAR(2,7) = 0D0
        C_X_STAR(3,7) = -7046897949D0/250000000000D0*SQRT(7D0)

!       Column 8 of C_X_STAR
        C_X_STAR(1,8) = 0D0
        C_X_STAR(2,8) = 0D0
        C_X_STAR(3,8) = 0D0

        C_TAU_STAR = 0.4D0
        C_P = 6
    END SELECT

  END SUBROUTINE INTERP_TABLEAU
      
    
  SUBROUTINE INTERP_WEIGHTS(S_STAR,TAU,W,WP)
! This routine evaluates the weight polynomials, and optionally their 
! first derivatives, at the relative position TAU with a subinterval. The
! weight polynomials correspond to the S_STAR stages used to form the 
! continuous Runge-Kutta interpolant. The results are returned in W and WP. 
!----------------------------------------------------------------------------
!     CALLED BY: DEFECT_ESTIMATE,INTERP_EVAL,SOL_EVAL
!----------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: S_STAR
    DOUBLE PRECISION :: TAU

!   Output arguments:
    DOUBLE PRECISION :: W(S_STAR)
    DOUBLE PRECISION, OPTIONAL :: WP(S_STAR)
  
!   Global variables:
!   BVP_METHOD identifies which discrete and continuous Runge-Kutta methods are 
!     being used, and thus which polynomials should be evaluated.
!----------------------------------------------------------------------------
    SELECT CASE (BVP_METHOD)
      CASE (2)
!       2-stage, 2nd order CMIRK formula.
        W(1) = TAU*(1D0 - TAU/2D0)
        W(2) = TAU**2 /2D0  
        IF ( PRESENT(WP) ) THEN  ! Derivative polynomials
          WP(1) = 1D0 - TAU
          WP(2) = TAU
        END IF
      CASE (4)
!       Improved 4-stage, 4th order CMIRK formula.
        W(1) = -TAU*(3D0*TAU-4D0)* (5D0*TAU**2 - 6D0*TAU + 3D0)/12D0
        W(2) = TAU**2*(5D0*TAU**2 - 6D0*TAU + 2D0)/6D0
        W(3) = -2D0/3D0*TAU**2*(3D0*TAU - 2D0)*(5D0*TAU - 6D0)
        W(4) = 125D0/12D0*TAU**2*(TAU - 1D0)**2
        IF ( PRESENT(WP) ) THEN  ! Derivative polynomials
          WP(1) = -(TAU-1.D0)*(5.D0*TAU-2.D0)*(2.D0*TAU-1.D0)/2.D0
          WP(2) = TAU*(5D0*TAU - 2D0)*(2D0*TAU - 1D0)/3D0
          WP(3) = -8D0*TAU*(TAU - 1D0)*(5D0*TAU - 2D0)
          WP(4) = 125D0/6D0*TAU*(TAU - 1D0)*(2D0*TAU - 1D0)
        END IF
      CASE (6)
!       Improved 8-stage, 6th order CMIRK formula. 
        W(1) =                                                                   &
        -(12233+1450*SQRT(7D0))*(800086000*TAU**5+63579600*                      &
        SQRT(7D0)*TAU**4-2936650584D0*TAU**4+4235152620D0*TAU**3                 &
        -201404565*SQRT(7D0)*TAU**3+232506630D0*SQRT(7D0)*TAU**2                 &
        -3033109390D0*TAU**2+1116511695D0*TAU-116253315D0*SQRT(7D0)              &
        *TAU+22707000D0*SQRT(7D0)-191568780D0)*TAU/2112984835740D0

        W(2) = -(-10799D0+650D0*SQRT(7D0))*(24962000D0*TAU**4+                   &
        473200D0*SQRT(7.D0)*TAU**3-67024328D0*TAU**3-751855D0*                   &
        SQRT(7D0)*TAU**2+66629600D0*TAU**2-29507250D0*TAU+236210*                &
        SQRT(7D0)*TAU+5080365D0+50895D0*SQRT(7D0))*TAU**2/                       &
        29551834260D0

        W(3) = 7D0/1274940D0*(259D0+50D0*SQRT(7D0))*                             &
        (14000D0*TAU**4-48216D0*TAU**3+1200D0*SQRT(7D0)*TAU**3-                  &
        3555D0*SQRT(7D0)*TAU**2+62790D0*TAU**2+3610D0*SQRT(7D0)*                 &
        TAU-37450D0*TAU+9135D0-1305D0*SQRT(7D0))*TAU**2

        W(4) = 7D0/1274940D0*(259D0+50D0*SQRT(7D0))*(14000D0*TAU**4              &
        -48216D0*TAU**3+1200D0*SQRT(7D0)*TAU**3-3555D0*SQRT(7D0)*                &
        TAU**2+62790D0*TAU**2+3610D0*SQRT(7D0)*TAU-37450D0*TAU+                  &
        9135D0-1305D0*SQRT(7D0))*TAU**2

        W(5) = 16D0/2231145D0*(259D0+50D0*SQRT(7D0))*(14000D0*                   &
        TAU**4- 48216D0*TAU**3+1200D0*SQRT(7D0)*TAU**3-3555D0*                   &
        SQRT(7D0)*TAU**2+62790D0*TAU**2+3610D0*SQRT(7D0)*                        &
        TAU-37450D0*TAU+9135D0-1305D0*SQRT(7D0))*TAU**2

        W(6) = 4D0/1227278493D0*(740D0*SQRT(7D0)-6083D0)*                        &
        (1561000D0*TAU**2-2461284D0*TAU-109520D0*SQRT(7D0)*TAU+                  &
        979272D0+86913D0*SQRT(7D0))*(TAU-1D0)**2*TAU**2

        W(7) = -49D0/63747D0*SQRT(7D0)*(20000D0*TAU**2-20000D0*                  &
        TAU+3393D0)*(TAU-1D0)**2*TAU**2

        W(8) = -1250000000D0/889206903D0*(28D0*TAU**2-28D0*TAU+                  &
        9D0)*(TAU-1D0)**2*TAU**2

        IF ( PRESENT(WP) ) THEN  ! Derivative polynomials
          WP(1) = (1450D0*SQRT(7D0)+12233D0)*(14D0*TAU-7D0+                      &
          SQRT(7D0))*(TAU-1D0)*(-400043D0*TAU+75481D0+2083D0*                    &
          SQRT(7D0))*(100D0*TAU-87D0)*(2D0*TAU-1D0)/493029795006D0

          WP(2) = -(650D0*SQRT(7D0)-10799D0)*(14D0*TAU-7D0+                      &
          SQRT(7D0))*(37443D0*TAU-13762D0-2083D0*SQRT(7D0))*                     &
          (100D0*TAU-87D0)*(2D0*TAU-1D0)*TAU/20686283982D0

          WP(3) = 7D0/42498D0*(259D0+50D0*SQRT(7D0))*(14D0*TAU-                  &
          7D0+SQRT(7D0))*(TAU-1D0)*(100D0*TAU-87D0)*(2D0*TAU-1D0)*TAU
          
          WP(4) = 7D0/42498D0*(259D0+50D0*SQRT(7D0))*(14D0*TAU-                  &
          7D0+SQRT(7D0))*(TAU-1D0)*(100D0*TAU-87D0)*(2D0*TAU-1D0)*TAU

          WP(5) = 32D0/148743D0*(259D0+50D0*SQRT(7D0))*(14D0*TAU-                &
          7D0+SQRT(7D0))*(TAU-1D0)*(100D0*TAU-87D0)*(2D0*TAU-1D0)*TAU

          WP(6) = 4D0/1227278493D0*(740D0*SQRT(7D0)-6083D0)*                     &
          (14D0*TAU-7D0+SQRT(7D0))*(TAU-1D0)*(100D0*TAU-87D0)*                   &
          (6690D0*TAU-4085D0-869D0*SQRT(7D0))*TAU

          WP(7) = -98D0/21249D0*SQRT(7D0)*(TAU-1D0)*(100D0*TAU-                  &
          13D0)*(100D0*TAU-87D0)*(2D0*TAU-1D0)*TAU

          WP(8) = -1250000000D0/2074816107D0*(14D0*TAU-7D0+SQRT(7D0))*           &
          (TAU-1D0)*(14D0*TAU-7D0-SQRT(7D0))*(2D0*TAU-1D0)*TAU
        END IF

    END SELECT

  END SUBROUTINE INTERP_WEIGHTS
  

  SUBROUTINE INTERVAL(NSUB,MESH,T,I)
! Finds I such that MESH(I-1) =< T < MESH(I).  Returns 1 if T <= MESH(0) 
! and NSUB if MESH(NSUB) <= T. MESH contains the current mesh of NSUB
! subintervals.
!---------------------------------------------------------
!     CALLED BY: INTERP_EVAL, SOL_EVAL
!----------------------------------------------------------
!   Input arguments:
    INTEGER :: NSUB
    DOUBLE PRECISION :: MESH(0:NSUB),T

!   Output arguments:
    INTEGER :: I

!   Local variables:
    INTEGER :: NDX(1)
!----------------------------------------------------------
    IF ( T <= MESH(0) ) THEN
      I = 1
    ELSE IF ( MESH(NSUB) <= T ) THEN
      I = NSUB
    ELSE
      NDX = MINLOC( MESH,MASK=(MESH >= T) )
      I = NDX(1) - 1   ! Indices start at 0 in MESH.
    END IF

  END SUBROUTINE INTERVAL

  SUBROUTINE JACBLK(NEQN,HI,TI,YI,YIP1,K,LLI,RRI,FSUB,DFSUB)
! This routine computes the left and right blocks of the Jacobian of the 
! residual function associated with the subinterval I. They are returned 
! in the vectors LLI and RRI.  Each is of size NEQN**2, where NEQN is 
! the number of ODEs plus the number of unknown parameters. The size of the 
! subinterval I is HI and its left endpoint is TI. The solution at the left 
! end of the subinterval is stored in YI and the solution at the right end is
! stored in YIP1. This computation uses the previously computed Runge-Kutta 
! stages provided in the vector K. FSUB and DFSUB are user-supplied routines
! for the ODEs and their derivatives.
!----------------------------------------------------------------------------
!     CALLED BY: NEWMAT
!     CALLS: P_DF,PD_F
!----------------------------------------------------------------------------
!   Modified by Mark Adams, 1996, to avoid unnecessary matrix multiplies.
!----------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NEQN
    DOUBLE PRECISION :: HI,TI,YI(NEQN),YIP1(NEQN),K(MXS*NEQN)

!   Output arguments:
    DOUBLE PRECISION :: LLI(NEQN**2),RRI(NEQN**2)

!   User-supplied subroutines:
    EXTERNAL FSUB,DFSUB
    
!   Local variables:
    INTEGER :: I,R,J,I1,I2 ! Loop and array indicies.
    DOUBLE PRECISION :: SUM ! Accumulator.
    INTEGER :: ALPHA, BETA, SIGMA ! Array indicies.      
    DOUBLE PRECISION :: TR ! T argument for Rth DF/DY evaluation.
            
    LOGICAL :: DKDYI_ZERO(MXS) ! Records presence of zeros in partial
    LOGICAL :: DKDYIP1_ZERO(MXS) ! derivatives to avoid matrix multiplies.
            
    LOGICAL :: MATRIX_MULT_DKDYI ! Indicator of special form of Runge-Kutta
    LOGICAL :: MATRIX_MULT_DKDYIP1 ! stage that can avoid matrix multiplies.
            
    DOUBLE PRECISION :: YR(NEQN) ! Y argument for Rth DF/DY evaluation.
    DOUBLE PRECISION :: PD(NEQN,NEQN) ! Store partial derivative.
    DOUBLE PRECISION :: DFDY(NEQN**2) ! Reshaped PD.
    DOUBLE PRECISION :: UR(NEQN**2) ! Intermediate values arising in 
                                    ! computation of partial derivative.                                           
            
    DOUBLE PRECISION :: DKDYI(NEQN**2*MXS)   ! Partial derivatives wrt
    DOUBLE PRECISION :: DKDYIP1(NEQN**2*MXS) ! YI and YIP1 associated with
                                             ! the derivatives of the Runge-
                                             ! Kutta stages.

!   Global variables:
!   MXS is the number of maximum number of Runge-Kutta stages.
!   NEQNSQ equals NEQN**2.
!   C_S is the number of discrete Runge-Kutta stages. 
!   The vectors C_C, C_V, C_B, and the array C_X contain the coefficients 
!   that define the discrete Runge-Kutta formula. 
!------------------------------------------------------------------------------
!   Compute partial derivatives in pairs.  In Rth pass through the loop, form 
!   d(KR)/d(YI) and d(KR)/d(YI+1) and store them in the appropriate segments of 
!   DKDYI and DKDYIP1.
!
!   Set DKDYI_ZERO and DKDYIP1_ZERO to .TRUE. to indicate that all of the partial  
!   derivatives are currently zero.
    DKDYI_ZERO   = .TRUE.
    DKDYIP1_ZERO = .TRUE.
      
    DO R = 1,C_S
!     PART 1: COMPUTE DF/DY. 
!     Define abscissa.
      TR = TI + C_C(R)*HI
      
!     Construct Y argument.
      YR = 0D0
!     Add on weighted contributions from first R-1 stages; stage J is  
!     available in the Jth set of NEQN locations of K.
      DO J = 1,R-1
        I1 = (J-1)*C_S + R
        I2 = (J-1)*NEQN
!       If C_X(R,J) is nonzero ...
        IF ( ABS(C_X(R,J)) > 0D0 ) THEN
          YR = YR + C_X(R,J)*K(I2+1:I2+NEQN)
        END IF
      END DO
      YR = (1D0 - C_V(R))*YI + C_V(R)*YIP1 + HI*YR
 
!     Arguments for DFDY are ready; compute DFDY by calling P_DF if analytic
!     form is available or by calling PD_F to approximate the derivatives 
!     using finite differences.
      IF (HAVE_DFDY) THEN
        CALL P_DF(DFSUB,TR,YR,PD,NEQN)
      ELSE
        CALL PD_F(FSUB,TR,YR,PD,NEQN)
      END IF

!     Output from P_DF or PD_F is a 2D array; reshape into vector DFDY.
      DFDY = RESHAPE(PD,(/NEQNSQ/))


!     PART 2: COMPUTE D(KR)/D(YI) AND STORE RESULT IN DKDYI.
!     Set MATRIX_MULT_DKDYI to .TRUE. indicating that the partial derivative 
!     is not of the form (1-C_V(R))*dF/dY and that multiplying the matrix UR
!     by DFDY may be necessary.
      MATRIX_MULT_DKDYI = .TRUE.
      
!     Accumulate weighted linear combination of previous R-1 partial 
!     derivative matrices in UR. 
      UR = 0D0
      DO J = 1, R-1
        I1 = (J-1)*C_S + R
        I2 = (J-1)*NEQNSQ

!       Skip if C_X(I1) = 0 or the Jth partial derivative is zero. 
        IF ( ABS(C_X(R,J)) > 0D0  .AND. ( .NOT. DKDYI_ZERO(J) ) ) THEN
          DKDYI_ZERO(R) = .FALSE. ! UR nonzero and we will need to do a 
                                  ! matrix multiply later.
          UR = UR + C_X(R,J)*DKDYI(I2+1:I2+NEQNSQ)
        END IF
      END DO
           
!     If UR is nonzero, scale it by HI.
      IF ( .NOT. DKDYI_ZERO(R) ) THEN
        UR = HI*UR
      END IF

      ALPHA = (R-1)*NEQNSQ
      IF ( ABS(C_V(R) - 1D0) > 0D0 ) THEN
        IF ( DKDYI_ZERO(R) ) THEN
!         Multiplying UR by DFDY will be unnecessary later in the computation.
          MATRIX_MULT_DKDYI = .FALSE.
          DKDYI(ALPHA+1:ALPHA+NEQNSQ) = ( 1D0 - C_V(R) )*DFDY
        ELSE
          UR(1:NEQNSQ:NEQN+1) = 1D0 - C_V(R) + UR(1:NEQNSQ:NEQN+1)
        END IF

!       Partial derivative of stage R is nonzero.
        DKDYI_ZERO(R) = .FALSE.
      END IF

!     Multiply UR and DFDY. Since both are stored as vectors, we code
!     the multiply with loops, dot products, and array index offsets.
!
!     Skip the multiplication if UR is zero or if partial derivative is
!     just (1-C_V(R))*DFDY.
      IF ( MATRIX_MULT_DKDYI ) THEN
        IF ( .NOT. DKDYI_ZERO(R) ) THEN
          DO I =1, NEQN
            BETA = ALPHA + I
            DO J = 1,NEQN
              SIGMA = (J-1)*NEQN
              DKDYI(SIGMA+BETA) = DOT_PRODUCT(DFDY(I:(NEQN-1)*NEQN+I:NEQN),      &
                                              UR(SIGMA+1:SIGMA+NEQN) )
            END DO
          END DO
        ELSE
!         If UR is zero, then the Rth partial derivative is zero. Value of 
!         DKDYI_ZERO(R) will left as .TRUE. to avoid later matrix multiplies.
          DKDYI(ALPHA+1:ALPHA+NEQNSQ) = 0D0
        END IF ! .NOT. DKDYI_ZERO(R) 
      END IF ! MATRIX_MULT_DKDYI


!     PART 3: COMPUTE D(KR)/D(YI+1) AND STORE RESULT IN DKDYIP1.
!     Set MATRIX_MULT_DKDYIP1 to .TRUE. This indicates that the partial derivative 
!     is not of the form C_V(R)*dF/dY and that multiplying the matrix UR by
!     DFDY may be necessary.
      MATRIX_MULT_DKDYIP1 = .TRUE.

!     Accumulate weighted linear combination of previous R-1 partial 
!     derivative matrices in UR.                  
      UR = 0D0
      DO J = 1, R-1
        I1 = (J-1)*C_S + R
        I2 = (J-1)*NEQNSQ
!       Skip if C_X(I1) is zero or the Jth partial derivative is zero.    
        IF ( ABS(C_X(R,J)) > 0D0  .AND.( .NOT. DKDYIP1_ZERO(J)) ) THEN
          DKDYIP1_ZERO(R) = .FALSE. ! UR is nonzero, so we will need to do a  
                                    ! matrix multiply later.
          UR = UR + C_X(R,J)*DKDYIP1(I2+1:I2+NEQNSQ)
        END IF
      END DO

!     If UR is nonzero, scale it.
      IF ( .NOT. DKDYIP1_ZERO(R) ) THEN
        UR = HI*UR
      END IF

      IF ( ABS(C_V(R)) > 0D0 ) THEN
        IF ( DKDYIP1_ZERO(R) ) THEN
!         Multiplying UR by DFDY will be unnecessary later in the computation.
          MATRIX_MULT_DKDYIP1 = .FALSE.
          DKDYIP1(ALPHA+1:ALPHA+NEQNSQ) = C_V(R)*DFDY
        ELSE
          UR(1:NEQNSQ:NEQN+1) = C_V(R) + UR(1:NEQNSQ:NEQN+1)
        END IF

!       Partial derivative of stage R stage is nonzero.
        DKDYIP1_ZERO(R) = .FALSE.
      END IF

!     Multiply UR by DFDY using loops and dot products since matrices are stored 
!     as vectors. Skip multiplication if UR = 0 or if partial derivative is just 
!     C_V(R)*DFDY.
      IF ( MATRIX_MULT_DKDYIP1 ) THEN
        IF ( .NOT. DKDYIP1_ZERO(R) ) THEN
          DO I = 1,NEQN
            BETA = ALPHA + I
            DO J = 1,NEQN
              SIGMA = (J-1)*NEQN
              DKDYIP1(SIGMA+BETA) = DOT_PRODUCT(DFDY(I:(NEQN-1)*NEQN+I:NEQN),    &
                                                UR(SIGMA+1:SIGMA+NEQN) )
            END DO
          END DO
        ELSE
!         If UR = 0 then Rth partial derivative is 0. Value of DKDYIP1_ZERO(R)
!         will left as .TRUE. to avoid later matrix multiplies.
          DKDYIP1(ALPHA+1:ALPHA+NEQNSQ) = 0D0
        END IF ! .NOT. DKDYIP1_ZERO(R)
      END IF !  MATRIX_MULT_DKDYIP1
    END DO


!   PART 4: COMPUTATION OF THE BLOCKS, LLI AND RRI.
    LLI = 0D0
    RRI = 0D0
    DO R = 1,C_S
      I1 = (R-1)*NEQNSQ
      LLI = LLI - C_B(R)*DKDYI(I1+1:I1+NEQNSQ)
      RRI = RRI - C_B(R)*DKDYIP1(I1+1:I1+NEQNSQ)
    END DO
    LLI = HI*LLI
    RRI = HI*RRI
!   Add on -I and I, respectively
    LLI(1:NEQNSQ:NEQN+1) = -1D0 + LLI(1:NEQNSQ:NEQN+1)
    RRI(1:NEQNSQ:NEQN+1) =  1D0 + RRI(1:NEQNSQ:NEQN+1)

  END SUBROUTINE JACBLK


  SUBROUTINE MESH_SELECTOR(NEQN,NSUB,MESH_CURRENT,DEFECT,NSUB_STAR,MESH_NEW,INFO)
! The purpose of this routine is to determine a new mesh. The NSUB+1 points 
! describing the current mesh are provided in MESH_CURRENT. An estimate of the 
! scaled defect of each solution component on each subinterval of the current 
! mesh is input through DEFECT. NEQN is the number of differential equations 
! plus the number of unknown parameters. This routine determines NSUB_STAR, 
! the number of subintervals for the new mesh from accuracy considerations. 
! The placement of new mesh points approximately equidistributes a monitor 
! function based on the scaled defect. The new mesh is returned in MESH_NEW. 
! If NSUB_STAR becomes greater than MXNSUB, the routine fails and returns
! with INFO = -1. Success is indicated by INFO = 0.
!--------------------------------------------------------------------------------
!     CALLED BY: MIRKDC
!     CALLS : HALF_MESH,REDISTRIBUTE
!--------------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NEQN,NSUB
    DOUBLE PRECISION :: MESH_CURRENT(0:NSUB),DEFECT(NEQN*NSUB)

!   Output arguments:
    INTEGER :: NSUB_STAR,INFO
    DOUBLE PRECISION MESH_NEW(0:4*NSUB) ! New mesh will have at most 4 times as
                                        ! many subintervals as current mesh.
!   Local variables:
    INTEGER :: I ! Loop index.
    INTEGER :: OFFSET ! Array index.
    INTEGER :: NSUB_PRED ! Predicted number of subintervals in new mesh.
    DOUBLE PRECISION :: HI ! Size of subinterval I of current mesh.
    DOUBLE PRECISION :: R1 ! Maximum size of weighted mesh function values.                                  
    DOUBLE PRECISION :: R2 ! Sum of weighted mesh function values; it gives an 
                           ! estimate of predicted number of points for new mesh.
    DOUBLE PRECISION :: R3 ! Average of weighted mesh function values.
    DOUBLE PRECISION :: NORM ! Max norm of defect for subinterval I.
    DOUBLE PRECISION :: SAFETY_FACTOR = 1.3D0 ! Safety factor to overestimate 
                                              ! number of points in new mesh.
    DOUBLE PRECISION :: RHO = 1D0 ! Threshold to trigger a remeshing based on ratio
                                  ! of maximum mesh function value to average value.
                                  ! (RHO = 1D0 forces a redistribution everytime.)
    DOUBLE PRECISION :: UPPER_NEW_MESH = 4D0   ! Upper and lower bounds on ratio.
    DOUBLE PRECISION :: LOWER_NEW_MESH = 0.5D0 ! of NSUB_STAR to NSUB.
    DOUBLE PRECISION :: S_HAT(NSUB) ! Mesh function values, one component per 
                                    ! subinterval. 
!   Global variables:
!   PROFILE gives the level of output to be provided.
!   C_P is the order of the Runge-Kutta formula.  
!   MXNSUB is the maximum number of subintervals.
!   BVP_TOL is the user-provided tolerance applied to the defect.
!-----------------------------------------------------------------------------
    INFO = 0
    R1 = 0D0
    R2 = 0D0

    DO I = 1,NSUB
      HI = MESH_CURRENT(I) - MESH_CURRENT(I-1)
      OFFSET = (I-1)*NEQN

!     Compute norm of scaled defect for subinterval I.
      NORM = MAXVAL(ABS(DEFECT(OFFSET+1:OFFSET+NEQN)))

!     Compute corresponding mesh function value.
      S_HAT(I) = (NORM/BVP_TOL)**(1D0/(C_P+1)) / HI

!     Update maximum weighted mesh function value, if necessary.
      R1 = MAX(R1,S_HAT(I)*HI)

!     Accumulate Ith weighted mesh function value in R2.
      R2 = R2 + S_HAT(I)*HI
    END DO

!   Compute average of weighted mesh function values.
    R3 = R2/NSUB

!   Estimate number of subintervals that new mesh should have. 
    NSUB_PRED = SAFETY_FACTOR*R2 + 1

!   To avoid cycling, make sure new mesh has a few more points than previous one.
    IF ( ABS( (NSUB_PRED - NSUB)/(1D0*NSUB)) < 0.1D0) NSUB_PRED = 1.1D0 * NSUB

!   Compute new mesh of NSUB_PRED points if maximum weighted mesh function value
!   is more than RHO times the average value. Otherwise simply halve the mesh.
    IF (R1 <= RHO*R3) THEN
!     Halve current mesh.
      NSUB_STAR = 2*NSUB
      IF (NSUB_STAR > MXNSUB) THEN
!       New mesh would be too large.
        INFO = -1
        IF (PROFILE > 0) THEN
          PRINT *,'Number of subintervals needed for new mesh would ',           &
          'exceed current allowed maximum of ', MXNSUB,'.'
          PRINT *,' '
        END IF
      ELSE 
        IF (PROFILE > 0) THEN
          PRINT *,'Halve each subinterval of the current mesh and try again.'
          PRINT *,'The number of subintervals in the new mesh will be ', NSUB,'.'
          PRINT *,' '
        END IF
        CALL HALF_MESH(NSUB,MESH_NEW)
      END IF ! (NSUB_STAR > MXNSUB)  
    ELSE
!     Compute new mesh of NSUB_PRED points.  
      NSUB_star = NSUB_pred
      IF ( NSUB_STAR > UPPER_NEW_MESH * NSUB ) NSUB_STAR = UPPER_NEW_MESH * NSUB
      IF ( NSUB_STAR < LOWER_NEW_MESH * NSUB ) NSUB_STAR = LOWER_NEW_MESH * NSUB
      IF (NSUB_STAR > MXNSUB) THEN
!       New mesh would be too large.
        INFO = -1
        IF (PROFILE > 0) THEN
          PRINT *,'Number of subintervals needed for new mesh would ',           &
          'exceed current allowed maximum of ', MXNSUB,'.'
          PRINT *,' '
        END IF
      ELSE 
        CALL REDISTRIBUTE(MESH_CURRENT,S_HAT,NSUB_STAR,MESH_NEW)          
      END IF ! (NSUB_STAR > MXNSUB)
    END IF ! (R1 <= RHO*R3)

  END SUBROUTINE MESH_SELECTOR
  

  SUBROUTINE NEWITER(NEQN,NSUB,MESH,Y,INFO,K_DISCRETE,FSUB,GSUB,DFSUB,DGSUB)
! This routine takes an approximate solution Y on a mesh, stored in MESH, of NSUB
! subintervals and computes a new approximation by applying a modified Newton's 
! method to the nonlinear system PHI(Y) = 0, where PHI(Y) is the residual 
! function. If the computation succeeds, the new approximation is returned in Y.
! As a by-product, the Runge-Kutta stage values are returned in the array 
! K_DISCRETE. The subroutines FSUB,GSUB,DFSUB, and DGSUB define the BVP. NEQN is
! the number of ODEs plus the number of unknown parameters. If the computation 
! succeeds, INFO = 0.  It can fail in three ways: If the Newton iteration takes 
! too many steps, INFO = 1. If the Newton matrix is singular, INFO = 2. If a 
! suitable damping factor cannot be found, INFO = 3. FSUB,GSUB,DFSUB,DGSUB are 
! user-supplied subroutines for the ODEs and boundary conditions and their 
! derivatives.
!-------------------------------------------------------------------------------
!      CALLED BY: MIRKDC
!      CALLS: FIXED_JACOB,DAMPED_NEWT,RESID
!-------------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NEQN, NSUB
    DOUBLE PRECISION :: MESH(0:NSUB),Y(NEQN*(NSUB+1))

!   Output arguments:
    INTEGER :: INFO
    DOUBLE PRECISION :: K_DISCRETE(MXS*NEQN*NSUB)

!   User-supplied subroutines:
    EXTERNAL FSUB,GSUB,DFSUB,DGSUB

!   Local variables:
    DOUBLE PRECISION :: PHI(NEQN*(NSUB+1)) ! Residual function.
    DOUBLE PRECISION :: DELTA(NEQN*(NSUB+1)) ! Newton correction.
    DOUBLE PRECISION :: TOP(NEQN**2) ! Derivatives of left boundary conditions.
    DOUBLE PRECISION :: BOT(NEQN**2) ! Derivatives of right boundary conditions.
    DOUBLE PRECISION  :: BLOCKS(2*NEQN**2 * NSUB) ! Derivatives of PHI(Y).
!   Together TOP, BLOCKS, and BOT contain an almost block diagonal matrix 
!   that is the Newton matrix. After a damped Newton step, TOP, BOT, 
!   and BLOCKS contain the factored form of this matrix.
    INTEGER :: PIVOT(NEQN*(NSUB+1)) ! Pivoting information from the 
                                    ! factorization of the Newton matrix.
    LOGICAL :: CONVRG ! Monitors convergence of Newton iteration.
    INTEGER :: COUNT ! Counts number of damped Newton iterations.
    DOUBLE PRECISION :: DELTA_0_NORM ! Norm of Newton correction vector.
    DOUBLE PRECISION :: LAMBDA ! Damping factor for damped Newton step.
    DOUBLE PRECISION :: G ! Natural criterion function value.
    LOGICAL :: FIXED_JACOBIAN ! .TRUE. if next step is a fixed Jacobian step; 
                              ! .FALSE. if next step is a damped Newton step.
    INTEGER, PARAMETER :: MAXITER=40 ! Maximum number of damped Newton steps.

!   Global variables:
!   MXS is the number of maximum number of Runge-Kutta stages.
!   PROFILE specifies the level of output.
!-----------------------------------------------------------------------------
!   Initialization for modified Newton iteration loop.
    COUNT = 0
    CONVRG = .FALSE.
    INFO = 0

!   The first iteration step must be a damped Newton step.
    FIXED_JACOBIAN = .FALSE.

!   When a damped Newton step is used, it is assumed that the current iterate 
!   Y and residual function PHI evaluated at Y are available. For the first 
!   step, Y is input, but PHI must be initialized through a call to the RESID 
!   routine.
    
    CALL RESID(NEQN,NSUB,MESH,Y,PHI,K_DISCRETE,FSUB,GSUB)
    
!
!   In predicting a new damping factor, certain quantities from the previous 
!   iteration are used. LAMBDA = 0 signals that they are not available for
!   the first step. 
    LAMBDA = 0D0

!   Hybrid damped Newton/Fixed Jacobian iteration. Report failure with INFO > 0.
    DO WHILE ( (.NOT. CONVRG) .AND.  (INFO <= 0) )
      IF (FIXED_JACOBIAN) THEN  
!       Don't count Fixed Jacobian steps--they are low cost.
  
        IF (PROFILE > 1) THEN
          PRINT *,'FIXED JACOBIAN STEP:' 
          PRINT *,' '
        END IF
        
        CALL FIXED_JACOB(NEQN,NSUB,MESH,Y,DELTA,G,PHI,TOP,BOT,BLOCKS,PIVOT,      &
                         FIXED_JACOBIAN,LAMBDA,CONVRG,INFO,K_DISCRETE,FSUB,GSUB)
      ELSE
        COUNT = COUNT + 1
       
        IF (PROFILE > 1) THEN
          PRINT *,'DAMPED NEWTON STEP, FULL NEWTON ITERATION COUNT =',COUNT,'.'
          PRINT *,' '
        END IF
        
        CALL DAMPED_NEWT(NEQN,NSUB,MESH,Y,LAMBDA,PHI,TOP,BOT,BLOCKS,PIVOT,       &
                         FIXED_JACOBIAN,CONVRG,DELTA,DELTA_0_NORM,G,INFO,        &
                         K_DISCRETE,FSUB,GSUB,DFSUB,DGSUB)
      END IF 

      IF ((.NOT. CONVRG) .AND. (INFO <= 0)) THEN
        IF ( COUNT > MAXITER) INFO = 1
      END IF
     
    END DO 

    IF (PROFILE > 0) THEN
      IF (COUNT > 1) THEN 
        WRITE(*,"(1X,A5,I3,A13)") 'AFTER',COUNT,' NEWTON STEPS,'
        PRINT *,' '
      ELSE
        WRITE(*,"(1X,A5,I3,A12)") 'AFTER',COUNT,' NEWTON STEP,'
        PRINT *,' '
      END IF
    END IF

  END SUBROUTINE NEWITER


  SUBROUTINE NEWMAT(LEFTBC,RIGHTBC,NEQN,NSUB,MESH,Y,TOP,BLOCKS,BOT,              &
                    K_DISCRETE,FSUB,GSUB,DFSUB,DGSUB)
! This routine evaluates the Newton matrix at the current iterate Y and 
! returns it in TOP, BOT, and BLOCKS; TOP contains the partial derivatives 
! associated with the LEFTBC boundary conditions; BOT contains the corresponding 
! information for the RIGHTBC boundary conditions; BLOCKS contains partial
! derivatives of the residual function with respect to the unknowns, y_i, 
! i=1,..., NSUB, where NSUB is the number of subintervals of the current MESH. 
! TOP stores a LEFTBC by NEQN matrix in column form; BOT stores a RIGHTBC by 
! NEQN matrix, also in column form. BLOCKS stores pairs of matrices L_i and R_i, 
! partial derivatives associated with subinterval i. Each of R_i and L_i is an 
! NEQN by NEQN matrix stored in column form. NEQN is the number of ODEs plus 
! the number of unknown parameters. This routine accesses the discrete 
! Runge-Kutta stages through the array K_DISCRETE which was assigned values in 
! a previous call to RESID with the same Y. FSUB,GSUB,DFSUB,DGSUB are 
! user-supplied subroutines for the ODEs and boundary conditions and their 
! derivatives.
!----------------------------------------------------------------------------
!      CALLED BY: DAMPED_NEWT
!      CALLS: JACBLK
!-----------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: LEFTBC,RIGHTBC,NEQN,NSUB
    DOUBLE PRECISION :: MESH(0:NSUB),Y(NEQN*(NSUB+1))
    DOUBLE PRECISION :: K_DISCRETE(MXS*NEQN*NSUB)

!   Output arguments:
    DOUBLE PRECISION ::  TOP(LEFTBC*NEQN),BLOCKS(2*NEQN**2*(NSUB+1)),            &
                         BOT(RIGHTBC*NEQN)

!   User-supplied subroutines:
    EXTERNAL FSUB,GSUB,DFSUB,DGSUB

!   Local variables:
    INTEGER :: I ! Loop index.       
    INTEGER :: I1,I2,I3,I4,N1 ! Array indicies.
    INTEGER :: TOPLOC,BOTLOC,WORKLOC ! Array indicies.
    DOUBLE PRECISION :: PD(NEQN,NEQN) ! Local storage of partial derivatives
                                      ! of boundary conditions.
    DOUBLE PRECISION :: WORK(NEQN**2) ! Local storage for reshaping of PD for
                                      ! transfer into TOP and BOT.
!   Global variables:
!   MXS is the number of maximum number of Runge-Kutta stages.
!   C_S is the number of discrete stages.
!   NEQNSQ equals NEQN**2.
!   HAVE_DBCDY specifies whether or not the user routine DBCDY was provided.
!-----------------------------------------------------------------------------
    I1 = C_S*NEQN 
    N1= NSUB*NEQN

    DO I = 1,NSUB
!     Compute blocks L_I and R_I through call to JACBLK. 
      I2 = (I-1)*NEQN ! Index offset for Y.
      I3 = (I-1)*I1 ! Index offset for K_DISCRETE.
      I4 = 2*(I-1)*NEQNSQ ! Index offset for BLOCKS.

!     y_I is stored starting at Y((I-1)*NEQN+1); y_{I+1} is stored starting at
!     Y(I*NEQN+1). The C_S stages associated with subinterval I are in the 
!     Ith set of C_S*NEQN locations of K_DISCRETE. L_I and R_I are returned in 
!     the Ith set of 2*NEQNS**2 locations of BLOCKS.
      CALL JACBLK(NEQN,MESH(I)-MESH(I-1),MESH(I-1),Y(I2+1),Y(I2+NEQN+1),         &
                  K_DISCRETE(I3+1),BLOCKS(I4+1),BLOCKS(I4+NEQNSQ+1),FSUB,DFSUB)
    END DO

!   If DGSUB is present, call P_DBC to get partial derivatives; otherwise
!   call PD_BC to compute the partial derivatives using finite differences.
!   Y(1:NEQN) is y_0, Y(N1+1:N1+NEQN) is y_NSUB. 
    IF (HAVE_DBCDY) THEN
      CALL P_DBC(DGSUB,Y(1:NEQN), Y(N1+1:N1+NEQN),PD(1:LEFTBC,1:NEQN),           &
                 PD(LEFTBC+1:NEQN,1:NEQN),LEFTBC,RIGHTBC,NEQN)               
    ELSE
      CALL PD_BC(GSUB,Y(1:NODE),Y(N1+1:N1+NEQN),PD(1:LEFTBC,1:NEQN),             &
                 PD(LEFTBC+1:NEQN,1:NEQN),LEFTBC,RIGHTBC,NEQN)
    END IF

    WORK = RESHAPE(PD,(/NEQNSQ/))

!   Since WORK contains the elements of the boundary condition matrix stored 
!   in a columnwise fashion, select the first LEFTBC elements from each set of 
!   NEQN elements of WORK and copy them into the appropriate locations of TOP. 
!   Similarly, copy the bottom RIGHTBC elements of each set of NEQN elements 
!   of WORK into the appropriate locations of BOT.

    DO I = 1,NEQN
      TOPLOC = (I-1)*LEFTBC
      WORKLOC = (I-1)*NEQN
      BOTLOC = (I-1)*RIGHTBC
      TOP(TOPLOC+1:TOPLOC+LEFTBC) = WORK(WORKLOC+1:WORKLOC+LEFTBC)
      BOT(BOTLOC+1:BOTLOC+RIGHTBC) = WORK(WORKLOC+LEFTBC+1:WORKLOC+NEQN)
    END DO

  END SUBROUTINE NEWMAT

  
  SUBROUTINE REDISTRIBUTE(MESH_CURRENT,S_HAT,NSUB_STAR,MESH_NEW)
! This routine is called with a mesh, MESH_CURRENT, of NSUB subintervals and a 
! vector of values, S_HAT, of a mesh function that is associated with a measure  
! of the relative defect. (More details about the mesh function are found in the
! routine MESH_SELECTOR.) It returns a new mesh, MESH_NEW, of NSUB_STAR 
! subintervals which equidistributes the integral of the mesh function. 
! The idea is to select new mesh points so that the area under the mesh function 
! is equally distributed over the subintervals of the new mesh. The quantity ZETA 
! is the total area under the mesh function divided by NSUB_STAR. The integral  
! over each subinterval of the new mesh is to equal ZETA. The area for each  
! subinterval of the old mesh is computed and added to a sum, proceeding from  
! left to right. Whenever the sum exceeds ZETA, the fraction of the subinterval  
! needed to equal ZETA is added. This fraction is used to place the new mesh 
! point. The procedure is repeated until every subinterval of the current mesh is
! processed. Chapter 9 of the text by Ascher, Mattheij, and Russell provides 
! further details.
!------------------------------------------------------------------------------
!   CALLED BY: MESH_SELECTOR
!------------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NSUB_STAR
    DOUBLE PRECISION :: MESH_CURRENT(0:NSUB),S_HAT(NSUB)

!   Output arguments:
    DOUBLE PRECISION :: MESH_NEW(0:NSUB_STAR)

!   Local variables:
    INTEGER :: I ! Index for tracking new mesh points.
    INTEGER :: K ! Index of subinterval of current mesh whose corresponding
                 ! mesh function integral is being added to integral for the
                 ! subinterval associated with the next new mesh point.
    DOUBLE PRECISION :: ZETA ! See above.
    DOUBLE PRECISION :: T ! Marker to keep track of how much of mesh function
                          ! integral has been processed.
    DOUBLE PRECISION :: INTEGRAL ! Keeps track of accumuation of integral to be
                                 ! associated with subinterval corresponding to
                                 ! next new mesh point. 0 <= INTEGRAL <= ZETA.
    DOUBLE PRECISION :: NEXT_PIECE ! Portion of mesh function integral processed
                                   ! during each interation of while loop.
!   Global variables:
!   NSUB is the number of subintervals in the current mesh.
!-------------------------------------------------------------------------------
!   Compute ZETA, the total area under the mesh function divided by NSUB_STAR.
    ZETA = SUM(S_HAT*(MESH_CURRENT(1:NSUB)-MESH_CURRENT(0:NSUB-1)))/NSUB_STAR
!
!   The first entry in both meshes is the left end point of the interval.    
    MESH_NEW(0) = MESH_CURRENT(0)

!   Set indicies, markers, and accumulators.
    K = 1
    I = 0
    T = MESH_CURRENT(0)
    INTEGRAL = 0D0

    DO WHILE (K <= NSUB)
      
      NEXT_PIECE = S_HAT(K)*(MESH_CURRENT(K) - T)
      
      IF ( (INTEGRAL + NEXT_PIECE) > ZETA) THEN
!       Define new mesh point so that the integral for this new subinterval is 
!       equal to ZETA.
        MESH_NEW(I+1) = (ZETA - INTEGRAL)/S_HAT(K) + T

!       Setup for next new meshpoint.
        T = MESH_NEW(I+1)
        INTEGRAL = 0D0        
        I = I + 1
      ELSE
!       Contribute all of remaining part of integral from subinterval K of 
!       current mesh to integral for the subinterval associated with next 
!       new mesh point.
        INTEGRAL = INTEGRAL + NEXT_PIECE
        
!       Set T and K to refer to next subinterval of current mesh.
        T = MESH_CURRENT(K)
        K = K + 1
      END IF

    END DO

!   The last entry in both meshes is the right end point of the interval.
    MESH_NEW(NSUB_STAR) = MESH_CURRENT(NSUB)
                                                                              
  END SUBROUTINE REDISTRIBUTE


  SUBROUTINE RESID(NEQN,NSUB,MESH,Y,PHI,K_DISCRETE,FSUB,GSUB)
! This subroutine evaluates the residual function PHI at the current iterate Y.  
! The current mesh of NSUB subintervals is stored in MESH. Evaluation of PHI 
! is done in two phases. The components associated with the interior of the 
! interval are evaluated through calls to the subroutine SUBCOM.  The components 
! associated with the boundary conditions are then obtained by by calls to 
! P_GSUB. A by-product of the call to SUBCOM is the return of the S Runge-Kutta 
! stages that arise as intermediate quantities. Each stage is a vector of length 
! NEQN where NEQN equals the number of ODEs plus the number of unknown 
! parameters. The stages are stored in order in the K_DISCRETE vector. FSUB,GSUB 
! are user-supplied subroutines for the ODEs and boundary conditions.
!-------------------------------------------------------------------------------
!     CALLED BY: NEWITER,CRITERION,FIXED_JACOB
!     CALLS: SUBCOM,P_GSUB
!-------------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NEQN,NSUB
    DOUBLE PRECISION :: MESH(0:NSUB),Y(NEQN*(NSUB+1))

!   Output arguments:
    DOUBLE PRECISION :: PHI(NEQN*(NSUB+1)),K_DISCRETE(MXS*NEQN*NSUB)

!   User-supplied subroutines:
    EXTERNAL FSUB,GSUB

!   Local variables:
    INTEGER :: I ! Loop index.
    INTEGER :: I1,N1 ! Array indicies.
    DOUBLE PRECISION :: H,T ! Subinterval length and left endpoint of 
                            ! subinterval.
!   Global variables:
!   MXS is the maximum number of Runge-Kutta stages.
!   C_S is the number of discrete stages.
!   LEFTBC is the number of left boundary conditions.
!-------------------------------------------------------------------------------
!   For each subinterval of the mesh,
    DO I = 1,NSUB
      H = MESH(I) - MESH(I-1)
      T = MESH(I-1)
      I1 = (I - 1)*NEQN
      
!     The current approximate solutions at the end points of subinterval I are
!     required for the computation of the residual.  The solution at the left
!     end is stored starting at Y((I-1)*NEQN+1) and the solution at the right
!     end is stored starting at Y(I*NEQN+1).  In the call to SUBCOM where the
!     residual is evaluated, PHI(LEFTBC+(I-1)*NEQN+1) is the segment of
!     PHI where this component of the residual is to be stored.  Similarly, the
!     stages of the discrete formula are returned in a segment of K_DISCRETE.
      
      CALL SUBCOM(NEQN,C_S,H,T,Y(I1+1:I1+NEQN),                                  &
                  Y(I1+NEQN+1:I1+2*NEQN),PHI(LEFTBC+I1+1:LEFTBC+I1+NEQN),        &
                  K_DISCRETE(I1*C_S+1:I1*C_S+NEQN*C_S),FSUB)

    END DO

!   The approximate solution at the left end point is stored starting at
!   Y(1) and the approximate solution at the right end point is stored
!   starting at Y(NSUB*NEQN+1). Results are returned in the first LEFTBC
!   locations of PHI and the last NEQN-LEFTBC locations of PHI.
    N1 = NSUB*NEQN
    CALL P_GSUB(GSUB,NEQN,LEFTBC,Y(1:NEQN),Y(N1+1:N1+NEQN),                      &
                PHI(1:LEFTBC),PHI(N1+LEFTBC+1:N1+NEQN))

  END SUBROUTINE RESID


  SUBROUTINE RK_TABLEAU
! This routine sets up the coefficients for the discrete mono-implicit 
! Runge-Kutta formula according to the value of the global variable BVP_METHOD. 
! The coefficients are assigned to global variables and arrays as indicated 
! below.
!-------------------------------------------------------------------------------
!     CALLED BY: BVP_SOLVER
!     CALLS: CHECK_STAT
!-------------------------------------------------------------------------------
!   Local variable: 
    INTEGER IER ! Error flag for allocation call.
       
!   Global variables:
!   BVP_METHOD identifies the Runge-Kutta method to be set up.
!   C_S is the number of Runge-Kutta stages. 
!   C_C, C_V, C_X, and C_B are arrays containing the Runge-Kutta coefficients.
!------------------------------------------------------------------------------------
    ALLOCATE(C_C(C_S),C_V(C_S),C_B(C_S),C_X(C_S,C_S),STAT=IER)
    CALL CHECK_STAT(IER)

    SELECT CASE (BVP_METHOD)
      CASE (2)
!       Define 2-stage, 2nd order, stage order 1 MIRK formula.
!       C_S has been set to 2.

        C_C(1) = 0D0
        C_C(2) = 1D0

        C_V(1) = 0D0
        C_V(2) = 1D0

        C_B(1) = 1D0/2D0
        C_B(2) = 1D0/2D0

!       Column 1 of C_X
        C_X(1,1) = 0D0
        C_X(2,1) = 0D0

!       Column 2 of C_X
        C_X(1,2) = 0D0
        C_X(2,2) = 0D0
      CASE (4)
!       Define 3-stage, 4th order, stage order 3 MIRK formula.
!       C_S has been set to 3.
      
        C_C(1) = 0D0
        C_C(2) = 1D0
        C_C(3) = 1D0/2D0

        C_V(1) = 0D0
        C_V(2) = 1D0
        C_V(3) = 1D0/2D0

        C_B(1) = 1D0/6D0
        C_B(2) = 1D0/6D0
        C_B(3) = 2D0/3D0

!       Column 1 of C_X
        C_X(1,1) = 0D0
        C_X(2,1) = 0D0
        C_X(3,1) = 1D0/8D0

!       Column 2 of C_X
        C_X(1,2) = 0D0
        C_X(2,2) = 0D0
        C_X(3,2) = -1D0/8D0

!       Column 3 of C_X
        C_X(1,3) = 0D0
        C_X(2,3) = 0D0
        C_X(3,3) = 0D0
      CASE (6)
!       Define improved 5-stage, 6th order, stage order 3 MIRK formula.
!       C_S has been set to 5.

        C_C(1) = 0D0
        C_C(2) = 1D0
        C_C(3) = 1D0/2D0 - SQRT(21D0)/14D0
        C_C(4) = 1D0/2D0 + SQRT(21D0)/14D0
        C_C(5) = 1D0/2D0

        C_V(1) = 0D0
        C_V(2) = 1D0
        C_V(3) = 1D0/2D0 - 9D0*SQRT(21D0)/98D0
        C_V(4) = 1D0/2D0 + 9D0*SQRT(21D0)/98D0
        C_V(5) = 1D0/2D0

        C_B(1) = 1D0/20D0
        C_B(2) = 1D0/20D0
        C_B(3) = 49D0/180D0
        C_B(4) = 49D0/180D0
        C_B(5) = 16D0/45D0

!       Column 1 of C_X          
        C_X(1,1) = 0D0
        C_X(2,1) = 0D0
        C_X(3,1) = 1D0/14D0 + SQRT(21D0)/98D0
        C_X(4,1) = 1D0/14D0 - SQRT(21D0)/98D0
        C_X(5,1) = -5D0/128D0

!       Column 2 of C_X
        C_X(1,2) = 0D0
        C_X(2,2) = 0D0
        C_X(3,2) = -1D0/14D0 + SQRT(21D0)/98D0
        C_X(4,2) = -1D0/14D0 - SQRT(21D0)/98D0
        C_X(5,2) = 5D0/128D0
!
!       Column 3 of C_X
        C_X(1,3) = 0D0
        C_X(2,3) = 0D0
        C_X(3,3) = 0D0
        C_X(4,3) = 0D0
        C_X(5,3) = 7D0*SQRT(21D0)/128D0
!
!       Column 4 of C_X
        C_X(1,4) = 0D0
        C_X(2,4) = 0D0
        C_X(3,4) = 0D0
        C_X(4,4) = 0D0
        C_X(5,4) = -7D0*SQRT(21D0)/128D0
         
!       Column 5 of C_X
        C_X(1,5) = 0D0
        C_X(2,5) = 0D0
        C_X(3,5) = 0D0
        C_X(4,5) = 0D0
        C_X(5,5) = 0D0
    END SELECT
   
  END SUBROUTINE RK_TABLEAU

    
  SUBROUTINE SOL_EVAL(NODE,NEQN,IWORK,WORK,T,Z,Z_PRIME)
! This routine evaluates the interpolant at T.  The value is returned in Z and
! optionally the value of the first derivative is returned in Z_PRIME. NEQN is 
! NODE, the number of ODEs, plus the number of unknown parameters. IWORK and WORK 
! contain information needed by the interpolant. The IWORK array contains
! S, the number of discrete Runge-Kutta stages, S_STAR, the total number of 
! Runge-Kutta stages, and BVP_METHOD, which specifies which Runge-Kutta method is
! used. The WORK array contains K_DISCRETE, the discrete stages; K_INTERP, the 
! continuous stages; MESH, the mesh points; and Y, the current solution 
! approximation at the mesh points.
!--------------------------------------------------------------------------------
!     CALLED BY: BVP_EVAL
!     CALLS: INTERVAL,INTERP_WEIGHTS,SUM_STAGES
!--------------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NODE,NEQN,IWORK(:)
    DOUBLE PRECISION :: WORK(:),T

!   Output arguments:
    DOUBLE PRECISION :: Z(NODE)
    DOUBLE PRECISION, OPTIONAL :: Z_PRIME(NODE)

!   Local variables:
    INTEGER :: S,S_STAR ! Number of discrete and continous Runge-Kutta stages.
    INTEGER :: I ! Index to Ith subinterval.
    INTEGER :: I_K_DISCRETE,I_K_INTERP,I_MESH,I_Y ! Work array indices.
    DOUBLE PRECISION :: HI ! Size of subinterval I.
    DOUBLE PRECISION :: TAU ! Relative position of T in subinterval I.
    DOUBLE PRECISION :: WEIGHTS(MXS) ! Weight polynomials for interpolant
                                     ! evaluated at TAU.
    DOUBLE PRECISION :: WEIGHTS_PRIME(MXS) ! Derivatives of weight polynomials
                                           ! evaluated at TAU.
    LOGICAL :: HAVE_Z_PRIME ! Records whether or not Z_PRIME is present
    DOUBLE PRECISION :: Z_LONG(NEQN) ! Contains approximate solution at T in its
                                     ! first NODE locations and approximation to
                                     ! unknown parameters in its last NPAR 
                                     ! locations.
    DOUBLE PRECISION :: Z_PRIME_LONG(NEQN) ! Contains approximation to first 
                                           ! derivative at T in its first NODE 
                                           ! locations and zero in its last 
                                           ! NPAR locations.  
!   Global variables:
!   MXS is the maximum number of Runge-Kutta stages.
!   BVP_METHOD identifies the Runge-Kutta methods.
!   NSUB is the number of subintervals of the mesh.
!--------------------------------------------------------------------------------
    HAVE_Z_PRIME = PRESENT(Z_PRIME)

!   Extract and assign values from IWORK.
    S = IWORK(1)
    S_STAR = IWORK(2)
    BVP_METHOD = IWORK(3) ! Assign value to global variable BVP_METHOD for use 
                          ! inside call to INTERP_WEIGHTS.

!   Set array indicies for convenient access to WORK.
    I_K_DISCRETE = 1
    I_K_INTERP = S*NEQN*NSUB+1
    I_MESH = S_STAR*NEQN*NSUB+1
    I_Y = S_STAR*NEQN*NSUB+(NSUB+1)+1

!   The INTERVAL subroutine returns the index of an interval containing T,
!   if T is within the problem domain. Extrapolation is used for T outside 
!   [A,B] = [WORK(I_MESH), WORK(I_MESH+NSUB)].
    CALL INTERVAL(NSUB,WORK(I_MESH:I_MESH+NSUB),T,I)
                                                                        
!   Compute TAU, the relative location of T in subinterval I
!   for use in the evaluation of the weight polynomials.
    HI = WORK(I_MESH+I) - WORK(I_MESH+(I-1))
    TAU = (T - WORK(I_MESH+(I-1)))/HI

!   Evaluate WEIGHTS and optionally WEIGHTS_PRIME at TAU.
    IF (HAVE_Z_PRIME) THEN
      CALL INTERP_WEIGHTS(S_STAR,TAU,WEIGHTS,WEIGHTS_PRIME)
    ELSE
      CALL INTERP_WEIGHTS(S_STAR,TAU,WEIGHTS)
    END IF

!   Evaluate interpolant at TAU, and optionally its derivative, by taking 
!   weighted sums of Runge-Kutta stages through a call to SUM_STAGES. Results are
!   returned in Z_LONG and Z_PRIME_LONG. The first NODE locations of Z_LONG 
!   contain the interpolant value; the remaining locations contain the 
!   approximation to the unknown parameters. The first NODE locations of 
!   Z_PRIME_LONG contain the derivative approximation; the remaining locations 
!   are zero.
    IF (HAVE_Z_PRIME) THEN
      CALL SUM_STAGES(NEQN,HI,WORK(I_Y+(I-1)*NEQN:I_Y+I*NEQN),S,                 &
                      WORK(I_K_DISCRETE+(I-1)*S*NEQN:I_K_DISCRETE+I*S*NEQN),     &
                      S_STAR,WORK(I_K_INTERP+(I-1)*(S_STAR-S)*NEQN:              &
                      I_K_INTERP+I*(S_STAR-S)*NEQN),WEIGHTS,                     &
                      Z_LONG,WEIGHTS_PRIME,Z_PRIME_LONG) 
    ELSE
      CALL SUM_STAGES(NEQN,HI,WORK(I_Y+(I-1)*NEQN:I_Y+I*NEQN),S,                 &
                      WORK(I_K_DISCRETE+(I-1)*S*NEQN:I_K_DISCRETE+I*S*NEQN),     &
                      S_STAR,WORK(I_K_INTERP+(I-1)*(S_STAR-S)*NEQN:              &
                      I_K_INTERP+I*(S_STAR-S)*NEQN),WEIGHTS,Z_LONG) 
    END IF

!   Copy first NODE locations of Z_LONG to Z; if Z_PRIME is present, copy 
!   first NODE locations of Z_PRIME_LONG to Z_PRIME.
    Z = Z_LONG(1:NODE)
    IF (HAVE_Z_PRIME) THEN
      Z_PRIME = Z_PRIME_LONG(1:NODE)
    END IF
    
  END SUBROUTINE SOL_EVAL
  

  SUBROUTINE SUBCOM(NEQN,C_S,HI,TI,YI,YIP1,PHII,K,FSUB)
! This subroutine evaluates the residual PHII associated with subinterval I.
! The basic formula is PHII = YIP1 - YI - HI* SUM BR*KR.  The evaluation
! takes place for the current iterate Y. However, component I of PHI depends 
! only on the solution information directly associated with subinterval I,
! namely YI and YIP1. The length of this subinterval is HI and the left hand
! end point is TI. During the computation of PHII, we compute C_S intermediate 
! vectors, KR, R = 1,..,C_S, called the stages of the formula. Each stage is a
! vector of length NEQN; NEQN equals the number of ODEs plus the number of 
! unknown parameters. The stages are returned in the vector K, one after the 
! other. The stages represent evaluations of the ODEs using the FSUB subroutine. 
!------------------------------------------------------------------------------
!   CALLED BY: RESID
!   CALLS: P_FSUB
!------------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NEQN,C_S
    DOUBLE PRECISION :: HI,TI,YI(NEQN),YIP1(NEQN)

!   Output arguments:
    DOUBLE PRECISION :: PHII(NEQN),K(C_S*NEQN)

!   User-supplied subroutines:
    EXTERNAL FSUB
           
!   Local variables:
    INTEGER :: R ! Loop index.
    DOUBLE PRECISION :: TR,Y_ARG(NEQN) ! Local storage for T and Y arguments
                                       ! in call to compute stage values.
    DOUBLE PRECISION :: KMATRIX(NEQN,C_S) ! Local representation of K as a
                                          ! 2D array. Reshaped to give the 1D
                                          ! vector K at end of computation.
!   Global variables:
!   The calculation uses the C_S-stage MIRK formula defined by  
!   C_S, C_C, C_V, C_B, and C_X.
!------------------------------------------------------------------------------
    KMATRIX = 0D0
    DO R = 1,C_S
      Y_ARG = (1D0 - C_V(R))*YI + C_V(R)*YIP1 + HI*MATMUL(KMATRIX,C_X(R,:))
      TR = TI + C_C(R)*HI
      CALL P_FSUB(FSUB,NEQN,TR,Y_ARG,KMATRIX(:,R))
    END DO
    
    PHII = YIP1 - YI - HI*MATMUL(KMATRIX,C_B)

    K = RESHAPE(KMATRIX,(/NEQN*C_S/))
               
  END SUBROUTINE SUBCOM


  SUBROUTINE SUM_STAGES(NEQN,HI,YIM1,S,KI_DISCRETE,S_STAR,KI_INTERP,WEIGHTS,Z,   &
                        WEIGHTS_PRIME,Z_PRIME)
! Given the stages and corresponding weights, this routine computes the value of 
! the continuous extension, Z, and (optionally) its first derivative, Z_PRIME,
! as follows:  Z = YIM1 + HI* SUM WEIGHTS(R)*KI(R)
!              Z_PRIME = SUM WEIGHTS_PRIME(R)*KI(R)
! Here the KI(R)'s are the stages.  The first S stages, contained in KI_DISCRETE,
! come from the discrete formula; the remaining (S_STAR-S) stages, contained in
! KI_INTERP, are necessary for the interpolant.  Each stage is a vector of
! length NEQN; NEQN equals the number of ODEs plus the number of unknown 
! parameters.  HI is the size of subinterval I and YIM1 is the discrete solution
! at the left end point of the subinterval.  WEIGHTS contains the values of the
! weight polynomials of the interpolant and WEIGHTS_PRIME contains the values of
! the first derivative of the weight polynomials.
!---------------------------------------------------------------------------
!   CALLED BY: DEFECT_EXTIMATE,INTERP_EVAL,SOL_EVAL
!---------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NEQN,S,S_STAR
    DOUBLE PRECISION :: HI,YIM1(NEQN),KI_DISCRETE(S*NEQN),                      &
                        KI_INTERP((S_STAR-S)*NEQN)
    DOUBLE PRECISION :: WEIGHTS(S_STAR)
    DOUBLE PRECISION, OPTIONAL :: WEIGHTS_PRIME(S_STAR)

!   Output arguments:
    DOUBLE PRECISION :: Z(NEQN)
    DOUBLE PRECISION, OPTIONAL :: Z_PRIME(NEQN)

!   Local variables:
    INTEGER :: J ! Loop index.
!---------------------------------------------------------------------------
!   Compute Z.
    Z = 0D0
    DO J = 1,S
      Z = Z + WEIGHTS(J)*KI_DISCRETE((J-1)*NEQN+1:J*NEQN+1)
    END DO
    DO J = 1,(S_STAR-S)
      Z = Z + WEIGHTS(S+J)*KI_INTERP((J-1)*NEQN+1:J*NEQN+1)
    END DO
    Z = YIM1 + HI*Z

!   If present, compute Z_PRIME.
    IF (PRESENT(Z_PRIME)) THEN
      Z_PRIME = 0D0
      DO J = 1,S
        Z_PRIME = Z_PRIME + WEIGHTS_PRIME(J)*KI_DISCRETE((J-1)*NEQN+1:J*NEQN+1)
      END DO
      DO J = 1,(S_STAR-S)
        Z_PRIME = Z_PRIME + WEIGHTS_PRIME(S+J)*KI_INTERP((J-1)*NEQN+1:J*NEQN+1)
      END DO
    END IF    
  
  END SUBROUTINE SUM_STAGES


!================beginning of generic BVP_INIT==========================

  FUNCTION GUESS_1(NODE,LEFTBC,X,Y,P,MAX_NUM_SUBINTERVALS) RESULT(SOL)
! Expand constant guess for the solution. Load vector guess, Y, 
! for the solution into the columns of the Y field of the solution 
! structure, SOL, which has one column for each point of the initial
! mesh, X. The vector Y gives NODE initial guesses for the solution 
! components and these will be stored in the first NODE locations of
! each column of Y. LEFTBC is the number of left boundary conditions. 
! We also load the values NPTS (the number of meshpoints) and NODE 
! into SOL and initialize the INFO field to zero. P if present gives
! initial guesses for the unknown parmeters and these are copied into
! the corresponding field of SOL. We also set SOL%MXNSUB, the maximum number 
! of subintervals, depending on the MAX_NUM_SUBINTERVALS argument. If X 
! contains only the end points, we set a default mesh of 10 equally 
! spaced points.
!-------------------------------------------------------------------------------
!     CALLED BY: User's program
!     CALLS: CHECK_STAT,BVP_LINSPACE
!--------------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NODE,LEFTBC
    DOUBLE PRECISION, DIMENSION(:) :: X,Y(NODE)
    DOUBLE PRECISION, DIMENSION(:), OPTIONAL :: P 
    INTEGER, OPTIONAL :: MAX_NUM_SUBINTERVALS

!   Function output:
    TYPE(BVP_SOL) :: SOL  
    
!   Local variables:
    INTEGER :: NPTS,IER ! Number of meshpoints and memory allocation error flag.   
    LOGICAL :: DEFAULT_MESH ! Flag to signal default mesh.
!----------------------------------------------------------------------------------
!   Set NPTS and check endpoints.
    NPTS = SIZE(X)
    IF (X(1) >= X(NPTS)) THEN
      PRINT *,'The interval [A,B] must have A < B.'
      STOP
    END IF

!  If NPTS = 2, use default of 10 equally spaced points.      
   DEFAULT_MESH = (NPTS == 2)
   IF (DEFAULT_MESH) NPTS = 10
  
!  Create storage for allocatable fields of SOL called X and Y, 
!  where initial mesh and initial solution guess and initial 
!  parameter guess, if present, will be stored.
   ALLOCATE(SOL%X(NPTS),SOL%Y(NODE,NPTS),STAT=IER)
   CALL CHECK_STAT(IER)

!  Load X and Y fields; copy Y vector into each column of SOL%Y.
   IF (DEFAULT_MESH) THEN
     SOL%X = BVP_LINSPACE(X(1),X(2),NPTS)
   ELSE
     SOL%X = X
   END IF
   SOL%Y = SPREAD(Y,DIM=2,NCOPIES=NPTS)

!  Load NPAR and P fields if P present.
   IF (PRESENT(P)) THEN
     ALLOCATE(SOL%PARAMETERS(SIZE(P)),STAT=IER)
     CALL CHECK_STAT(IER)
     SOL%PARAMETERS = P
     SOL%NPAR = SIZE(P)
   ELSE
     SOL%NPAR = 0
   END IF
   
!  Load remaining fields.
   SOL%NPTS = NPTS
   SOL%NODE = NODE
   SOL%LEFTBC = LEFTBC
   SOL%INFO = 0

!  Set SOL%MXNSUB.
    IF (PRESENT(MAX_NUM_SUBINTERVALS)) THEN
      SOL%MXNSUB = MAX_NUM_SUBINTERVALS
    ELSE
      SOL%MXNSUB = 3000
    ENDIF
    
  END FUNCTION GUESS_1
  

  FUNCTION GUESS_2(NODE,LEFTBC,X,Y,P,MAX_NUM_SUBINTERVALS) RESULT(SOL)
! Load a matrix guess, Y, for the solution - one column for each of 
! the points of the initial mesh, X, into the columns of SOL%Y. Each
! column of Y gives NODE initial values for the solution components 
! which are stored in the first NODE locations of each column of
! SOL%Y. LEFTBC is the number of left boundary conditions. We also 
! load the values NPTS (the number of meshpoints) and NODE into SOL 
! and initialize the INFO field to zero. P if present gives initial 
! guesses for the unknown parmeters and these are copied into the 
! corresponding field of SOL. We also set SOL%MXNSUB, the maximum number 
! of subintervals, depending on the MAX_NUM_SUBINTERVALS argument. If X 
! contains only the end points, use a default mesh of 10 equally spaced 
! points. 
!-----------------------------------------------------------------------
!     CALLED BY: User's program
!     CALLS: CHECK_STAT,BVP_LINSPACE
!-----------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NODE,LEFTBC
    DOUBLE PRECISION, DIMENSION(:) :: X
    DOUBLE PRECISION, DIMENSION(NODE,SIZE(X)) :: Y
    DOUBLE PRECISION, DIMENSION(:), OPTIONAL :: P
    INTEGER, OPTIONAL :: MAX_NUM_SUBINTERVALS

!   Function output:
    TYPE(BVP_SOL) :: SOL 

!   Local variables:
    INTEGER :: NPTS,IER ! Number of meshpoints and error flag for allocation 
                        ! call.   
    LOGICAL :: DEFAULT_MESH ! Flag to signal default mesh. 
!----------------------------------------------------------------------------
!   Set NPTS and check endpoints.
    NPTS = SIZE(X)
    IF (X(1) >= X(NPTS)) THEN
      PRINT *,'The interval [A,B] must have A < B.'
      STOP
    END IF

!   If NPTS = 2, use default of 10 equally spaced points.      
    DEFAULT_MESH = (NPTS == 2)
    IF (DEFAULT_MESH) NPTS = 10    
        
!   Create storage for allocatable fields of SOL called X and Y, 
!   where initial mesh and initial solution guess will be stored.
    ALLOCATE(SOL%X(NPTS),SOL%Y(NODE,NPTS),STAT=IER)
    CALL CHECK_STAT(IER)
    
!   Load X and Y fields; copy Y array into SOL%Y.
    IF (DEFAULT_MESH) THEN
      SOL%X = BVP_LINSPACE(X(1),X(2),NPTS)
    ELSE
      SOL%X = X
    END IF
    SOL%Y = Y

!   Load NPAR and P fields if P present.
    IF (PRESENT(P)) THEN
    ALLOCATE(SOL%PARAMETERS(SIZE(P)),STAT=IER)
    CALL CHECK_STAT(IER)
    SOL%PARAMETERS = P
    SOL%NPAR = SIZE(P)
    ELSE
    SOL%NPAR = 0
    END IF

!   Load remaining fields.
    SOL%NPTS = NPTS
    SOL%NODE = NODE
    SOL%LEFTBC = LEFTBC
    SOL%INFO = 0

!   Set SOL%MXNSUB.
    IF (PRESENT(MAX_NUM_SUBINTERVALS)) THEN
      SOL%MXNSUB = MAX_NUM_SUBINTERVALS
    ELSE
      SOL%MXNSUB = 3000
    ENDIF
    
  END FUNCTION GUESS_2
  

  FUNCTION GUESS_3(NODE,LEFTBC,X,FCN,P,MAX_NUM_SUBINTERVALS) RESULT(SOL)
! Evaluate a subroutine FCN for the initial solution guess at the 
! points of the initial mesh, X, and load into the columns of the 
! Y field of the solution structure, SOL. The call to FCN will 
! provide NODE initial values for the solution components and these
! will be stored in the first NODE locations of each column of Y. 
! NODE is the number of ODEs. LEFTBC is the number of left boundary 
! conditions. We also load the values NPTS (the number of meshpoints) 
! and NODE into SOL and initialize the INFO field to zero. P if present 
! gives initial guesses for the unknown parmeters and these are copied into
! the corresponding field of SOL. We also set SOL%MXNSUB, the maximum number 
! of subintervals, depending on the MAX_NUM_SUBINTERVALS argument. If X 
! contains only the end points, use a default mesh of 10 equally spaced 
! points.
!----------------------------------------------------------------------------
!     CALLED BY: User's program
!     CALLS: CHECK_STAT,BVP_LINSPACE
!----------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: NODE,LEFTBC
    DOUBLE PRECISION, DIMENSION(:) :: X
    EXTERNAL FCN
    DOUBLE PRECISION, DIMENSION(:), OPTIONAL :: P
    INTEGER, OPTIONAL :: MAX_NUM_SUBINTERVALS

!   Function output:
    TYPE(BVP_SOL) :: SOL 

!   Local variables:
    INTEGER :: NPTS,IER,K ! Number of meshpoints, error flag for allocation 
                          ! call, loop index.   
    LOGICAL :: DEFAULT_MESH ! Flag to signal default mesh.   
!----------------------------------------------------------------------------
!   Set NPTS and check endpoints.
    NPTS = SIZE(X)
    IF (X(1) >= X(NPTS)) THEN
      PRINT *,'The interval [A,B] must have A < B.'
      STOP
    END IF

!   If NPTS = 2, use default of 10 equally spaced points.      
    DEFAULT_MESH = (NPTS == 2)
    IF (DEFAULT_MESH) NPTS = 10    
    
!   Create storage for allocatable fields of SOL called X and Y, 
!   where initial mesh and initial solution guess and initial 
!   parameter guess, if present, will be stored.
    ALLOCATE(SOL%X(NPTS),SOL%Y(NODE,NPTS),STAT=IER)
    CALL CHECK_STAT(IER)

!   Load X and Y fields; copy vector from evaluation of function FCN into 
!   first NODE locations of each column of SOL%Y. 
    IF (DEFAULT_MESH) THEN
      SOL%X = BVP_LINSPACE(X(1),X(2),NPTS)
    ELSE
      SOL%X = X
    END IF
    DO K = 1,NPTS
    CALL FCN(SOL%X(K),SOL%Y(1:NODE,K))
    END DO

!   Load NPAR and P fields if P present.
    IF (PRESENT(P)) THEN
      ALLOCATE(SOL%PARAMETERS(SIZE(P)),STAT=IER)
      CALL CHECK_STAT(IER)
      SOL%PARAMETERS = P
      SOL%NPAR = SIZE(P)
    ELSE
      SOL%NPAR = 0
    END IF
    
!   Load remaining fields.
    SOL%NPTS = NPTS
    SOL%NODE = NODE
    SOL%LEFTBC = LEFTBC
    SOL%INFO = 0

!   Set SOL%MXNSUB.
    IF (PRESENT(MAX_NUM_SUBINTERVALS)) THEN
      SOL%MXNSUB = MAX_NUM_SUBINTERVALS
    ELSE
      SOL%MXNSUB = 3000
    ENDIF

  END FUNCTION GUESS_3
  
!================end of generic BVP_INIT===========================

    
  SUBROUTINE P_DF(DFSUB,T,Y,PD,NEQN)
! Evaluate analytical partial derivative d FSUB/dY at (T,Y) or d FSUB/dY and
! d FSUB/dP at (T,Y,P) if FSUB depends on unknown parameters. (P is available in
! Y(NODE+1:NEQN)). If the problem has a singular term, handle this separately. 
! DFSUB is the user function which provides the partial derivatives. The partial 
! derivatives are returned by this routine in the 2D array PD, which is of size 
! NEQN by NEQN, where NEQN equals the number of ODEs plus the number of unknown 
! parameters.
!------------------------------------------------------------------------------
!     CALLED BY: JACBLK
!     CALLS: DFSUB
!------------------------------------------------------------------------------
!   Input arguments:
    EXTERNAL DFSUB  
    INTEGER :: NEQN  
    DOUBLE PRECISION :: T,Y(NEQN)
    
!   Output argument:    
    DOUBLE PRECISION :: PD(NEQN,NEQN)

!   Local variable:
    DOUBLE PRECISION :: A

!   Global variables:
!   NPAR is the number of unknown paramters.
!   NODE is the number of ODEs.
!   SINGULAR indicates whether or not the BVP has a singular term.
!   FSING and MSING are quantities associated witht the singular term of
!     the BVP, if it is present.
!   G_X is the array of mesh points.
!------------------------------------------------------------------------------
    PD = 0D0

!   Call DFSUB with appropriate parameter list, based on whether or not 
!   there are unknown parameters. This loads the partial derivatives into PD.
    IF (NPAR == 0) THEN
      CALL DFSUB(T,Y,PD)
    ELSE
!     Call DFSUB with additional P argument (P is available in Y(NODE+1:NEQN)); 
!     it returns d FSUB/dY and d FSUB/dP in the last two arguments.
      CALL DFSUB(T,Y(1:NODE),Y(NODE+1:NEQN),&
                 PD(1:NODE,1:NODE),PD(1:NODE,NODE+1:NEQN))
    END IF
    
!   If the problem has a singular term, handle the
!   evaluation of the partial derivatives as a special case.
    IF (SINGULAR) THEN
!     Extract left end point A from the global array G_X.          
      A = G_X(0)
!     If T is the left endpoint.
      IF (T <= A) THEN
        PD(1:NODE,1:NODE) = MATMUL(FSING,PD(1:NODE,1:NODE))
      ELSE
        PD(1:NODE,1:NODE) = MSING/(T-A) + PD(1:NODE,1:NODE)
      END IF
    END IF            

  END SUBROUTINE P_DF


  SUBROUTINE PD_F(FSUB,T,Y,PD,NEQN)
! Compute divided difference approximation to d FSUB/dY at (T,Y) or d FSUB/dY and
! d FSUB/dP at (T,Y,P) if FSUB depends on unknown parameters. (P is available in 
! Y(NODE+1:NEQN)). FSUB is the user function which evaluates the right hand side 
! of the ODEs. The partial derivatives are returned by this routine in the 2D 
! array PD, which is of size NEQN by NEQN, where NEQN equals the number of ODEs 
! plus the number of unknown parameters.
!------------------------------------------------------------------------------
!     CALLED BY: JACBLK
!     CALLS: P_FSUB
!--------------------------------------------------------------------------------
!   Input arguments:
    EXTERNAL FSUB  
    INTEGER :: NEQN  
    DOUBLE PRECISION :: T,Y(NEQN)

!   Output argument: 
    DOUBLE PRECISION :: PD(NEQN,NEQN)

!   Local variables:
    INTEGER :: K ! Loop index.
    DOUBLE PRECISION :: SAVEYK ! Kth component of Y.
    DOUBLE PRECISION :: DELYK ! Denominator for divided difference.
    DOUBLE PRECISION :: FTY(NEQN) ! Value of FSUB at Y.
    DOUBLE PRECISION :: FDELY(NEQN) ! Value of FSUB at Y + DELYK.

!   Global variables:
!   SQRTU is the square root of unit roundoff.
!--------------------------------------------------------------------------------
!   Evaluate FSUB for input T and Y. P_FSUB calls FSUB with the appropriate 
!   argument list depending on whether or not FSUB depends on unknown parameters.   
    CALL P_FSUB(FSUB,NEQN,T,Y,FTY)

!   Compute Jacobian a column at a time.
    DO K = 1,NEQN
     
      SAVEYK = Y(K)
      
!     Determine appropriate perturbation of Y(K) and perturb Y(K) by this amount.
      IF ( ABS(SAVEYK) > 0D0 ) THEN
           DELYK = SQRTU*ABS(SAVEYK)
      ELSE
        DELYK = SQRTU
      END IF
      Y(K) = SAVEYK + DELYK

!     Evaluate FSUB for input T and perturbed Y.
      CALL P_FSUB(FSUB,NEQN,T,Y,FDELY)

!     Compute divided difference giving Kth column of Jacobian. 
      PD(:,K) = (FDELY - FTY) / DELYK
      
!     Restore Y(K) value.
      Y(K) = SAVEYK

    END DO

  END SUBROUTINE PD_F


  SUBROUTINE P_DBC(DGSUB,YA,YB,PDYA,PDYB,LEFTBC,RIGHTBC,NEQN)
! Evaluate analytical partial derivative d GSUB/dYA at (A,YA) and d GSUB/dYB 
! at (B,YB) or d GSUB/dYA and d GSUB/dP at (A,YA,P) and d GSUB/dYB and d GSUB/dP 
! at (B,YB,P) if GSUB depends on unknown parameters. (P is available in 
! YA(NODE+1:NEQN)). DGSUB is the user function which provides the partial 
! derivatives. The partial derivatives are returned by this routine in the 2D 
! arrays PDYA abd PDYB, which are of size LEFTBC by NEQN and RIGHTBC by NEQN, 
! where LEFTBC is the number of left boundary conditions, RIGHTBC is the number 
! of right boundary conditions, and NEQN equals the number of ODEs plus the 
! number of unknown parameters.
!------------------------------------------------------------------------------
!     CALLED BY: NEWMAT
!------------------------------------------------------------------------------
!   Input arguments:
    EXTERNAL DGSUB
    INTEGER :: LEFTBC,RIGHTBC,NEQN
    DOUBLE PRECISION :: YA(NEQN),YB(NEQN)

!   Output arguments:
    DOUBLE PRECISION :: PDYA(LEFTBC,NEQN),PDYB(RIGHTBC,NEQN)

!   Global variables:
!   NODE is the number of ODEs.
!   NPAR is the number of unknown parameters.
!-----------------------------------------------------------------------------
    PDYA = 0D0
    PDYB = 0D0

!   Call DGSUB with appropriate argument list depending on whether or not there
!   are unknown parameters. 
    IF (NPAR == 0) THEN
      CALL DGSUB(YA,YB,PDYA,PDYB)
    ELSE
!     The third and fourth arguments of DGSUB return d GSUB/dYA at (A,YA) and 
!     d GSUB/dYB at (B,YB); P (stored in YA(NODE+1:NEQN)) is input as the fifth 
!     argument and the sixth and seventh arguments of DGSUB return d GSUB/dP at 
!     (A,YA,P) and d GSUB/dP at (B,YB,P).
      CALL DGSUB(YA(1:NODE),YB(1:NODE),PDYA(:,1:NODE),&
                 PDYB(:,1:NODE),YA(NODE+1:NEQN),&
                 PDYA(:,NODE+1:NEQN),PDYB(:,NODE+1:NEQN))
    END IF

  END SUBROUTINE P_DBC
  

  SUBROUTINE PD_BC(GSUB,YA,YB,PDYA,PDYB,LEFTBC,RIGHTBC,NEQN)
! Compute divided difference approximation to d GSUB/dYA at (A,YA) and d GSUB/dYB 
! at (B,YB) or d GSUB/dYA and d GSUB/dP at (A,YA,P) and d GSUB/dYB and d GSUB/dP 
! at (B,YB,P) if GSUB depends on unknown parameters. (P is available in 
! YA(NODE+1:NEQN)). GSUB is the user function which evaluates the boundary 
! conditions. The partial derivatives are returned by this routine in the 2D 
! arrays  PDYA and PDYB. YA is of length LEFTBC and YB is of length RIGHTBC, so 
! PDYA is LEFTBC by NEQN and PDYB is RIGHTBC by NEQN, where LEFTBC is the number
! of left boundary conditions, RIGHTBC is the number of right boundary conditions, 
! and NEQN equals the number of ODEs plus the number of unknown parameters.
!------------------------------------------------------------------------------
!     CALLED BY: NEWMAT
!     CALLS: P_GSUB
!------------------------------------------------------------------------------
!   Input arguments:
    EXTERNAL GSUB
    INTEGER :: LEFTBC,RIGHTBC,NEQN
    DOUBLE PRECISION :: YA(NEQN),YB(NEQN)

!   Output arguments:
    DOUBLE PRECISION :: PDYA(LEFTBC,NEQN),PDYB(RIGHTBC,NEQN)
    
!   Local variables:
    INTEGER :: K ! Loop index.
    DOUBLE PRECISION :: SAVEYA,SAVEYB ! Kth components of YA and YB.
    DOUBLE PRECISION :: DELYA,DELYB ! Denominators for divided differences.
    DOUBLE PRECISION :: BCA(LEFTBC),BCB(RIGHTBC) ! Boundary conditions at YA and YB.
    DOUBLE PRECISION :: PBCA(LEFTBC),PBCB(RIGHTBC) ! Boundary conditions at 
                                                   ! YA + DELYA and YB + DELYB.

!   Global variables:
!   SQRTU is the square root of unit roundoff.
!------------------------------------------------------------------------------    
!   Evaluate GSUB for input YA and YB. P_GSUB calls GSUB with the appropriate 
!   argument list depending on whether or not GSUB depends on unknown parameters.   
    CALL P_GSUB(GSUB,NEQN,LEFTBC,YA,YB,BCA,BCB)

!   Compute Jacobian a column at a time.
    DO K = 1,NEQN   

      SAVEYA = YA(K)

!     Determine appropriate perturbation of YA(K) and perturb YA(K) by this amount.      
      IF ( ABS(SAVEYA) > 0D0 ) THEN
            DELYA = SQRTU*ABS(SAVEYA)
      ELSE
        DELYA = SQRTU
      END IF
      YA(K) = SAVEYA + DELYA

      SAVEYB = YB(K)

!     Determine appropriate perturbation of YB(K) and perturb YB(K) by this amount.      
      IF ( ABS(SAVEYB) > 0D0 ) THEN
            DELYB = SQRTU*ABS(SAVEYB)
      ELSE
        DELYB = SQRTU
      END IF
      YB(K) = SAVEYB + DELYB

!     Evaluate GSUB at perturbed YA and YB.      
      CALL P_GSUB(GSUB,NEQN,LEFTBC,YA,YB,PBCA,PBCB)

!     Compute divided difference giving the Kth columns of the Jacobians. 
      PDYA(:,K) = (PBCA - BCA) / DELYA
      PDYB(:,K) = (PBCB - BCB) / DELYB

!     Restore YA(K) and YB(K) values.
      YA(K) = SAVEYA
      YB(K) = SAVEYB

    END DO

  END SUBROUTINE PD_BC
      

  SUBROUTINE CHECK_STAT(IER)
! Check status of an attempted array allocation or deallocation by
! checking error flag IER. Print error condition and stop if 
! allocation/deallocation fails.
!---------------------------------------------------------------------
!   CALLED BY: BVP_SOLVER,MIRKDC,INTERP_TABLEAU,RK_TABLEAU,BVP_INIT,
!              BVP_GET,BVP_EXTEND
!---------------------------------------------------------------------
!   Input argument:
    INTEGER :: IER
!---------------------------------------------------------------------
      IF (IER /= 0) THEN
          PRINT *,'A storage allocation or deallocation error occurred.'
          STOP
      END IF
        
  END SUBROUTINE CHECK_STAT    
      
  
  FUNCTION BVP_LINSPACE(A,B,NPTS) RESULT(X)
! Generate a uniform mesh of NPTS points between A and B (inclusive);
! return mesh in array X, locations 1 to NPTS.
!---------------------------------------------------------------------
!   CALLED BY: User's program,BVP_INIT
!---------------------------------------------------------------------
!   Input arguments: 
    INTEGER :: NPTS
    DOUBLE PRECISION :: A,B
    
!   Function result:
    DOUBLE PRECISION ::  X(NPTS)

!   Local variables:
    INTEGER :: K ! Loop index.
    INTEGER :: MNPTS ! Local copy of NPTS
    DOUBLE PRECISION :: DELTA ! Subinterval size.
!----------------------------------------------------------------------
!   Mesh must have at least two points.
!   (The following line causes trouble with sample problems econ, lubrication,
!   and SWAVE, with some compilers unless we !   use a local copy of NPTS. 
!   Fix by Jason Boisvert/Gang Liang 2010.)
    MNPTS = MAX(NPTS,2)

    DELTA = (B - A) / (MNPTS - 1)
    X = (/ (A + DELTA*(K-1), K = 1,MNPTS) /)

  END FUNCTION BVP_LINSPACE
      

!================beginning of generic BVP_EVAL===========================

  SUBROUTINE EVAL_V(SOL,T,Z,DERIVATIVE)
! Evaluate solution (and optionaly derivative) at all the
! evaluation points contained in the vector T.
! Solutions returned in the columns of the 2D array Z; 
! derivatives in the columns of the 2D array DERIVATIVE,
! if this parameter is present. All info required for evaluation 
! of the solution and derivative is provided in SOL.
!------------------------------------------------------------------
!    CALLED BY: user's program 
!    CALLS: SOL_EVAL
!-----------------------------------------------------------------------
!  Input arguments: 
   TYPE(BVP_SOL) :: SOL
   DOUBLE PRECISION :: T(:)
   
!  Output arguments:
   DOUBLE PRECISION :: Z(SOL%NODE,SIZE(T))
   DOUBLE PRECISION, OPTIONAL :: DERIVATIVE(SOL%NODE,SIZE(T))

!  Local variables:
   INTEGER :: NEVALPTS ! Number of evaluation points, i.e. size of T.
   INTEGER :: K ! Loop index.

!  Global variables:
!  NSUB is the number of subintervals of the current mesh.
!-----------------------------------------------------------------------
    NEVALPTS = SIZE(T)
  
    IF (SOL%INFO /= 0) THEN
      WRITE(*,*) 'BVP_SOLVER failed, so the solution cannot be saved.'
      STOP
    END IF 
    
!   Set global variable NSUB; it is referenced inside SOL_EVAL.
    NSUB = SOL%NPTS - 1
    
    IF (PRESENT(DERIVATIVE)) THEN
      DO K = 1,NEVALPTS
        CALL SOL_EVAL(SOL%NODE,SOL%NODE+SOL%NPAR,SOL%IWORK,SOL%WORK,             &
                      T(K),Z(:,K),DERIVATIVE(:,K))
      END DO
    ELSE
      DO K = 1,NEVALPTS
        CALL SOL_EVAL(SOL%NODE,SOL%NODE+SOL%NPAR,SOL%IWORK,SOL%WORK,T(K),Z(:,K))
      END DO
    END IF        

  END SUBROUTINE EVAL_V

  
  SUBROUTINE EVAL_S(SOL,T,Z,DERIVATIVE)
! Evaluate solution (and optionaly derivative) at a single 
! point T. Solution returned in Z; derivative in DERIVATIVE, if 
! this parameter is present. All info required for evaluation
! of the solution and derivative in SOL.
!------------------------------------------------------------------
!    CALLED BY: user's program,BVP_EXTEND 
!    CALLS: SOL_EVAL
!-----------------------------------------------------------------------  
!  Input arguments: 
    TYPE(BVP_SOL) :: SOL
    DOUBLE PRECISION :: T

!  Output arguments:
   DOUBLE PRECISION :: Z(SOL%NODE)
   DOUBLE PRECISION, OPTIONAL :: DERIVATIVE(SOL%NODE)

!  Global variables:
!  NSUB is the number of subintervals of the current mesh.
!-----------------------------------------------------------------------
    IF (SOL%INFO /= 0) THEN
      WRITE(*,*) 'BVP_SOLVER failed, so the solution cannot be saved.'
      STOP
    END IF
    
!   Set global variable NSUB; it is referenced inside SOL_EVAL.
    NSUB = SOL%NPTS - 1

    IF (PRESENT(DERIVATIVE)) THEN
      CALL SOL_EVAL(SOL%NODE,SOL%NODE+SOL%NPAR,SOL%IWORK,SOL%WORK,T,Z,DERIVATIVE)
    ELSE
      CALL SOL_EVAL(SOL%NODE,SOL%NODE+SOL%NPAR,SOL%IWORK,SOL%WORK,T,Z)
    END IF        

  END SUBROUTINE EVAL_S


  SUBROUTINE EVAL_P(SOL,PARAMETERS)
! Return parameter values in PARAMETERS. All info required for evaluation
! of the parameters is provided in SOL.
!------------------------------------------------------------------
!    CALLED BY: user's program 
!-----------------------------------------------------------------------
!  Input arguments: 
   TYPE(BVP_SOL) :: SOL

!  Output arguments:    
   DOUBLE PRECISION :: PARAMETERS(SOL%NPAR)
!-----------------------------------------------------------------------
   IF (SOL%INFO /= 0) THEN
     WRITE(*,*) 'BVP_SOLVER failed, so the solution cannot be saved.'
     STOP
   END IF
    
   PARAMETERS = SOL%PARAMETERS
    
  END SUBROUTINE EVAL_P

!================end of generic BVP_EVAL===========================


  FUNCTION EYE(N) RESULT(IDENTITY)
! Returns a 2D array, IDENTITY, of size N by N containing identity matrix.
!-------------------------------------------------------------------------
!     CALLED BY: BVP_SOLVER
!-------------------------------------------------------------------------
!   Input argument:
    INTEGER :: N

!   Function result:
    DOUBLE PRECISION :: IDENTITY(N,N)

!   Local variables:
    INTEGER :: I,J ! Loop indexes.
!-------------------------------------------------------------------------
    IDENTITY = RESHAPE( (/ ((0D0,J=1,I-1),1D0,(0D0,J=I+1,N),                     &
                                   I=1,N) /), (/ N,N /) )        
  END FUNCTION EYE
 

  SUBROUTINE DUMMY_DF(T,Y,PD)
! Dummy subroutine used to fill in for missing DFSUB in call to MIRKDC 
! module. Has same argument list as DFSUB.
!-------------------------------------------------------------------------
!     CALLED BY: BVP_SOLVER
!-------------------------------------------------------------------------
!   Input arguments:
    DOUBLE PRECISION :: T,Y(:)
    
!   Output argument:    
    DOUBLE PRECISION :: PD(:,:)
!-------------------------------------------------------------------------
!   Dummy statements to avoid compiler warnings.
    T = 1
    Y(1) = 1
    PD(1,1) = 1

  END SUBROUTINE DUMMY_DF


  SUBROUTINE DUMMY_DBC(YA,YB,DYA,DYB)
! Dummy subroutine used to fill in for missing DGSUB in call to MIRKDC 
! module. Has same argument list as DGSUB.
!-------------------------------------------------------------------------
!     CALLED BY: BVP_SOLVER
!-------------------------------------------------------------------------
!   Input arguments:
    DOUBLE PRECISION :: YA(:),YB(:)

!   Output arguments:    
    DOUBLE PRECISION :: DYA(:,:),DYB(:,:)
!-------------------------------------------------------------------------
!   Dummy statements to avoid compiler warnings.
    YA(1) = 1
    YB(1) = 1
    DYA(1,1) = 1
    DYB(1,1) = 1

  END SUBROUTINE DUMMY_DBC


  SUBROUTINE P_FSUB(FSUB,NEQN,T,Y,FTY)
! Calls FSUB at (T,Y) or (T,Y,P) depending on whether or not FSUB depends 
! on unknown parameters, P. (P is available in Y(NODE+1:NEQN)). Result is 
! returned in FTY. NEQN is the number of ODEs plus the number of unknown 
! parameters. Also updates FSUB value to handle singular term, if it is present.
!------------------------------------------------------------------------
!     CALLED BY: PD_F,DEFECT_ESTIMATE,INTERP_SETUP,SUBCOM
!------------------------------------------------------------------------
!   Input arguments:
    EXTERNAL FSUB
    INTEGER :: NEQN
    DOUBLE PRECISION :: T,Y(NEQN)

!   Output argument:
    DOUBLE PRECISION :: FTY(NEQN)
    
!   Local variable:
    DOUBLE PRECISION :: A ! Left endpoint of problem interval.

!   Global variables:
!   NODE is the number of ODEs.
!   NPAR is the number of unknown parameters.
!   SINGULAR is TRUE if the BVP has a singular term; FALSE otherwise.
!   FSING,MSING are quantities associated with the treatment of the
!     singular term.
!   G_X contains the current mesh.
!------------------------------------------------------------------------
!   Call FSUB with appropriate parameter list, based on whether or not 
!   there are unknown parameters. This loads the function values into FTY.
    IF (NPAR == 0) THEN 
      CALL FSUB(T,Y,FTY)
    ELSE
!     Set FTY to zero since FSUB will only load first NODE locations.
!     Call FSUB with additional P argument (Y(NODE+1:NEQN)); function
!     values are loaded into first NODE locations of FTY.
      FTY = 0D0
      CALL FSUB(T,Y(1:NODE),Y(NODE+1:NEQN),FTY(1:NODE))
    END IF

!   If the problem has a singular term, handle the
!   evaluation of the function as a special case.
    IF (SINGULAR) THEN
!     Extract left end point A from the global array G_X.          
      A = G_X(0) 
!     If T is the left endpoint.
      IF (T <= A) THEN
        FTY(1:NODE) = MATMUL(FSING,FTY(1:NODE))
      ELSE
        FTY(1:NODE) = MATMUL(MSING,Y(1:NODE))/(T-A) + FTY(1:NODE)
      END IF
    END IF

  END SUBROUTINE P_FSUB


  SUBROUTINE P_GSUB(GSUB,NEQN,LEFTBC,YA,YB,BCA,BCB)
! Calls GSUB at (YA,YB) or (YA,YB,P) depending on whether or not GSUB depends 
! on unknown parameters, P. (P is available in Y(NODE+1:NEQN)). Results are 
! returned in BCA and BCB. NEQN is the number of ODEs plus the number of unknown 
! parameters. LEFTBC is the number of left boundary conditions.
!------------------------------------------------------------------------------------
!     CALLED BY: RESID,PD_BC
!------------------------------------------------------------------------------------ 
!   Input arguments:
    EXTERNAL GSUB
    INTEGER :: NEQN,LEFTBC
    DOUBLE PRECISION :: YA(NEQN),YB(NEQN)

!   Output arguments:
    DOUBLE PRECISION :: BCA(LEFTBC),BCB(NEQN-LEFTBC)

!   Global variables:
!   NPAR is the number of unknown parameters.
!   NODE is the number of ODEs.
!------------------------------------------------------------------------------------
!   Call DGSUB with appropriate argument list depending on whether or not there
!   are unknown parameters. 
    IF (NPAR == 0) THEN
      CALL GSUB(YA,YB,BCA,BCB)
    ELSE
!     The third argument is P=YA(NODE+1:NEQN).
      CALL GSUB(YA(1:NODE),YB(1:NODE),YA(NODE+1:NEQN),BCA,BCB)
    END IF

  END SUBROUTINE P_GSUB
 

  SUBROUTINE BVP_SAVE(UNUM,FNAME,SOL)
! BVP_SAVE saves the solution structure SOL in a file.  UNUM specifies 
! the unit to be OPENed.  The unit is CLOSEd after use. The string FNAME 
! specifies the name of the file. The structure can be retrieved from 
! this file with BVP_GET. Note that BVP_SAVE does not deallocate the
! array fields of SOL. This can be done afterwards with a call to
! BVP_TERMINATE.
!------------------------------------------------------------------------
!     CALLED BY: User's program.
!------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: UNUM
    CHARACTER(LEN=*) :: FNAME
    TYPE(BVP_SOL) :: SOL

!   Local variables:
    INTEGER :: LIWORK,LWORK ! Lengths of IWORK and WORK fields of SOL.
    INTEGER :: I,J ! Loop indexes.
!------------------------------------------------------------------------
    IF (SOL%INFO /= 0) THEN 
      WRITE(*,*) 'BVP_SOLVER failed, so the solution cannot be saved.'
      STOP
    END IF

    OPEN(UNIT=UNUM,FILE=FNAME,FORM="UNFORMATTED",                                &
         STATUS="REPLACE",POSITION="REWIND") 
    
    WRITE(UNUM) SOL%NODE,SOL%NPAR,SOL%LEFTBC,SOL%NPTS,SOL%INFO

!   Need the lengths of two arrays of pointers.
    LIWORK = SIZE(SOL%IWORK)
    LWORK = SIZE(SOL%WORK)
    WRITE(UNUM) LIWORK,LWORK
    
    DO I = 1,SOL%NPTS
      WRITE(UNUM) SOL%X(I)
    END DO
    DO I = 1,LIWORK
      WRITE(UNUM) SOL%IWORK(I)
    END DO
    DO I = 1,LWORK
      WRITE(UNUM) SOL%WORK(I)
    END DO
    DO I = 1,SOL%NODE
      DO J = 1,SOL%NPTS
        WRITE(UNUM) SOL%Y(I,J)
      END DO
    END DO
    IF (SOL%NPAR > 0) THEN
      DO I = 1,SOL%NPAR
        WRITE(UNUM) SOL%PARAMETERS(I)
      END DO
    END IF
    
    CLOSE(UNUM)

  END SUBROUTINE BVP_SAVE
  

  SUBROUTINE BVP_GET(UNUM,FNAME,SOL)
! BVP_GET forms the solution structure SOL from data stored by
! BVP_SAVE in the file FNAME.  UNUM specifies the unit to be 
! OPENed. The unit is CLOSEd after use.
!--------------------------------------------------------------------------
!     CALLED BY: User's program
!     CALLS: CHECK_STAT
!--------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: UNUM
    CHARACTER(LEN=*) :: FNAME

!   Output argument:    
    TYPE(BVP_SOL) :: SOL

!   Local variables:    
    INTEGER :: NODE,NPAR,NPTS ! Number of ODEs, number of unknown parameters,
                              ! Number of meshpoints.
    INTEGER :: LIWORK,LWORK ! Lengths of IWORK and WORK fields of SOL.
    INTEGER :: IER,I,J ! Error flag for array allocation, loop indexes.
!----------------------------------------------------------------------------
    OPEN(UNIT=UNUM,file=FNAME,FORM="UNFORMATTED",                                &
         STATUS="OLD",POSITION="REWIND") 
    READ(UNUM) SOL%NODE,SOL%NPAR,SOL%LEFTBC,SOL%NPTS,SOL%INFO
    NODE = SOL%NODE
    NPAR = SOL%NPAR    
    NPTS = SOL%NPTS
   
    READ(UNUM) LIWORK,LWORK
    
    IF (NPAR == 0) THEN
      ALLOCATE(SOL%X(NPTS),SOL%IWORK(LIWORK),SOL%WORK(LWORK),SOL%Y(NODE,NPTS),   &
              STAT=IER)
    ELSE
      ALLOCATE(SOL%X(NPTS),SOL%IWORK(LIWORK),SOL%WORK(LWORK),SOL%Y(NODE,NPTS),   &
              SOL%PARAMETERS(NPAR),STAT=IER)      
    END IF
    CALL CHECK_STAT(IER)
    
    DO I = 1,SOL%NPTS
      READ(UNUM) SOL%X(I)
    END DO
    DO I = 1,LIWORK
      READ(UNUM) SOL%IWORK(I)
    END DO
    DO I = 1,LWORK
      READ(UNUM) SOL%WORK(I)
    END DO
    DO I = 1,SOL%NODE
      DO J = 1,SOL%NPTS
        READ(UNUM) SOL%Y(I,J)
      END DO
    END DO
    IF (SOL%NPAR > 0) THEN
      DO I = 1,SOL%NPAR
        READ(UNUM) SOL%PARAMETERS(I)
      END DO
    END IF
   
    CLOSE(UNUM)

  END SUBROUTINE BVP_GET

  FUNCTION PINVSQ(A,N) RESULT(APLUS)
! Compute the pseudoinverse of an N by N matrix A. Result is returned in the
! N by N matrix APLUS. Uses LINPACK code DSVDC and its supporting programs in 
! BVP_LA.f. The tolerance is the same as in Matlab's PINV.
!-------------------------------------------------------------------------------
!     CALLED BY: BVP_SOLVER
!     CALLS: DSVDC
!-------------------------------------------------------------------------------
!   Input arguments:
    INTEGER :: N
    DOUBLE PRECISION :: A(N,N)

!   Function argument: 
    DOUBLE PRECISION :: APLUS(N,N)
   
!   Local variables:
    DOUBLE PRECISION :: COPYA(N,N) ! Saves copy of A.
    INTEGER :: INFO,I ! Info flag for DSVDC, Loop index.
    DOUBLE PRECISION :: S(N),U(N,N),V(N,N) ! S is a vector of singular values in 
                                           ! descending order of magnitude. U is 
                                           ! a matrix of left singular vectors. 
                                           ! V is a matrix of right singular 
                                           ! vectors.  
    DOUBLE PRECISION :: E(N),WORK(N) ! Work arrays required by DSVDC.
    DOUBLE PRECISION :: TOL ! Tolerance on the singular values.

!   Global variable:
!   UROUND is the machine unit round off.
    EXTERNAL DSVDC
!--------------------------------------------------------------------------------
!   Calculate the singular value decomposition with a call to DSVDC. DSVDC gives
!   S, U, V such that A = U * diag(S) * V'. Then A+ = V * diag(S)+ * U'.  
!   Here "+" indicates pseudoinverse. 

!   DSVDC destroys A which has the side effect of destroying it in 
!   the calling program, so use a local copy.
    COPYA = A
   
    CALL DSVDC(COPYA,N,N,N,S,E,U,N,V,N,WORK,11,INFO)
    
    IF (INFO /= 0) THEN
      PRINT *,'Computation of pseudoinverse failed.'
      STOP
    END IF
    
!   Form pseudoinverse with tolerance TOL on the singular values.
    TOL = N*S(1)*UROUND
    APLUS = 0D0
    DO I = 1,N
      IF (S(I) < TOL) EXIT
      APLUS(I,:) = U(:,I)/S(I)
    END DO
    APLUS = MATMUL(V,APLUS)  

 END FUNCTION PINVSQ


!================beginning of generic BVP_EXTEND==========================

  FUNCTION EXTEND_1(SOLIN,ANEW,BNEW,YANEW,YBNEW,P,MAX_NUM_SUBINTERVALS)          &
                                                        RESULT(SOLOUT)   
! This function allows the solution in a structure of type BVP_SOL to be
! extended to a solution on a bigger interval. A typical call has the form
!       SOLOUT = BVP_EXTEND(SOLIN,ANEW,BNEW,YANEW,YBNEW)
! Here SOLIN is the input solution structure and SOLOUT is the output 
! structure.  If the interval used to compute SOLIN was [A,B], the new 
! guess is for [ANEW,BNEW]. If ANEW < A, the mesh is extended to ANEW and 
! Y is extended to the value YANEW there. If ANEW >= A, the arguments ANEW
! and YANEW are ignored. The right end of the interval is treated in the 
! same way. If SOL%NPAR > 0, SOLIN%PARAMETERS is used in SOLOUT. To supply a 
! different guess P for unknown parameters, use the optional argument P.  
! To increase MXNSUB, the maximum number of subintervals allowed, use the 
! optional argument MAX_NUM_SUBINTERVALS.
!---------------------------------------------------------------------------
!     CALLED BY: User's program
!     CALLS: CHECK_STAT
!---------------------------------------------------------------------------
!   Input arguments:
    TYPE(BVP_SOL) :: SOLIN
    DOUBLE PRECISION :: ANEW,YANEW(:),BNEW,YBNEW(:)
    DOUBLE PRECISION, OPTIONAL :: P(:)
    INTEGER, OPTIONAL :: MAX_NUM_SUBINTERVALS

!   Function output:
    TYPE(BVP_SOL) :: SOLOUT

!   Local variables:
    INTEGER :: NODE,NPAR ! Number of ODEs, number of parameters.
    INTEGER :: LEFTBC ! Number of left boundary conditions.
    INTEGER :: NPTS ! Number of meshpoints. 
    INTEGER :: MXNSUB ! Maximum number of subintervals allowed.
    
    INTEGER :: IER ! Error flag for array allocation.
    INTEGER :: MODE ! Records type of extension to be performed.
    DOUBLE PRECISION :: A,B ! Problem interval endpoints.

!   Temporary storage of array info from SOLIN.
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: X,PARAMETERS
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: Y
!--------------------------------------------------------------------------
!   Assign scalar info from SOLIN structure to local variables.
    NODE = SOLIN%NODE
    NPAR = SOLIN%NPAR
    LEFTBC = SOLIN%LEFTBC
    NPTS = SOLIN%NPTS
    MXNSUB = SOLIN%MXNSUB

!   Allocate temporary storage for info from SOLIN array fields and assign
!   these fields to local arrays.
    ALLOCATE(X(NPTS),Y(NODE,NPTS),STAT=IER)
    CALL CHECK_STAT(IER)
    X = SOLIN%X
    Y = SOLIN%Y
    IF (NPAR > 0) THEN
      ALLOCATE(PARAMETERS(NPAR),STAT=IER)
      CALL CHECK_STAT(IER)
      PARAMETERS = SOLIN%PARAMETERS
    END IF

!   Relevant info from SOLIN is now unpacked. Note that IWORK and WORK 
!   arrays of SOLIN will no longer be valid after the mesh and solution 
!   info is extended. Deallocate all array fields of SOLIN.
    DEALLOCATE(SOLIN%X,SOLIN%Y,SOLIN%IWORK,SOLIN%WORK,STAT=IER)
    CALL CHECK_STAT(IER)
    IF (NPAR > 0) THEN
      DEALLOCATE(SOLIN%PARAMETERS,STAT=IER)
      CALL CHECK_STAT(IER)
    END IF


!   Determine required type of extension and set MODE and NPTS accordingly.
    A = X(1)
    B = X(NPTS)
    MODE = 0
    IF (ANEW < A) THEN
      NPTS = NPTS + 1
      MODE = 1
    END IF
    IF (B < BNEW) THEN
      NPTS = NPTS + 1      
      IF (MODE == 0) THEN
        MODE = 2
      ELSE
        MODE = 3
      END IF    
    END IF

!   Pack solution info into SOLOUT structure; allocate array fields as needed.
    SOLOUT%INFO = 0
    SOLOUT%NODE = NODE
    SOLOUT%NPAR = NPAR
    SOLOUT%LEFTBC = LEFTBC
    SOLOUT%NPTS = NPTS
    
    IF (PRESENT(MAX_NUM_SUBINTERVALS)) THEN
      SOLOUT%MXNSUB = MAX_NUM_SUBINTERVALS
    ELSE
      SOLOUT%MXNSUB = MXNSUB
    ENDIF

    IF (NPAR > 0) THEN
      ALLOCATE(SOLOUT%PARAMETERS(NPAR),STAT=IER)
      CALL CHECK_STAT(IER)
      SOLOUT%PARAMETERS = PARAMETERS
    END IF
    
!   If P is present, assign parameter values to SOLOUT field. Verify that
!   the number of parameters is not changed.
    IF (PRESENT(P)) THEN
      IF (NPAR /= SIZE(P)) THEN
        PRINT *,'Cannot change the number of unknown parameters.'
        STOP
      END IF
      SOLOUT%PARAMETERS = P
    END IF

!   Allocate arrays for mesh and solution in SOLOUT.
    ALLOCATE(SOLOUT%X(NPTS),SOLOUT%Y(NODE,NPTS),STAT=IER)
    CALL CHECK_STAT(IER)

!   Based on value of MODE, load mesh and solution fields appropriately.
    SELECT CASE(MODE)
      CASE (0) ! ANEW = A, B = BNEW
        SOLOUT%X = X
        SOLOUT%Y = Y
      CASE (1) ! ANEW < A, B = BNEW
        SOLOUT%X(1) = ANEW
        SOLOUT%X(2:NPTS) = X
        SOLOUT%Y(:,1) = YANEW
        SOLOUT%Y(:,2:NPTS) = Y
      CASE (2) ! ANEW = A, B < BNEW
        SOLOUT%X(1:NPTS-1) = X
        SOLOUT%X(NPTS) = BNEW
        SOLOUT%Y(:,1:NPTS-1) = Y
        SOLOUT%Y(:,NPTS) = YBNEW
      CASE (3) ! ANEW < A, B < BNEW
        SOLOUT%X(1) = ANEW
        SOLOUT%X(2:NPTS-1) = X
        SOLOUT%X(NPTS) = BNEW
        SOLOUT%Y(:,1) = YANEW
        SOLOUT%Y(:,2:NPTS-1) = Y
        SOLOUT%Y(:,NPTS) = YBNEW
    END SELECT

!   SOLOUT structure is now loaded. Deallocate temporary arrays.
    DEALLOCATE(X,Y,STAT=IER)
    CALL CHECK_STAT(IER)
    IF (NPAR > 0) THEN
      DEALLOCATE(PARAMETERS,STAT=IER)
      CALL CHECK_STAT(IER)
    END IF

  END FUNCTION EXTEND_1 
  

  FUNCTION EXTEND_2(SOLIN,ANEW,BNEW,ORDER,P,MAX_NUM_SUBINTERVALS)              &
                                                        RESULT(SOLOUT)   
! This function allows the solution in a structure of type BVP_SOL to be
! extended to a solution on a bigger interval. A typical call has the form
!       SOLOUT = BVP_EXTEND(SOLIN,ANEW,BNEW,ORDER)
! Here SOLIN is the input solution structure and SOLOUT is the output 
! structure.  If the interval used to compute SOLIN was [A,B], the new 
! guess is for [ANEW,BNEW]. If ANEW < A, the mesh is extended to ANEW and 
! Y is extended to a value YANEW there. If ANEW >= A, no action is taken
! at the left end.  The value YANEW is computed by extrapolation.  The 
! optional integer ORDER indicates the order of the polynomial used.  
! If ORDER = 1, linear extrapolation is done and otherwise, constant 
! extrapolation (ORDER = 0). The right end of the interval is treated in 
! the same way. If SOL%NPAR > 0, SOLIN%PARAMETERS is used in SOLOUT. To supply 
! a different guess P for unknown parameters, use the optional argument P.  
! To increase MXNSUB, the maximum number of subintervals allowed, use the 
! optional argument MAX_NUM_SUBINTERVALS.
!---------------------------------------------------------------------------
!     CALLED BY: User's program
!     CALLS: CHECK_STAT,BVP_EVAL
!---------------------------------------------------------------------------
!   Input arguments:
    TYPE(BVP_SOL) :: SOLIN
    DOUBLE PRECISION :: ANEW,BNEW
    INTEGER, OPTIONAL :: ORDER
    DOUBLE PRECISION, OPTIONAL :: P(:)
    INTEGER, OPTIONAL :: MAX_NUM_SUBINTERVALS

!   Function output:
    TYPE(BVP_SOL) :: SOLOUT

!   Local variables:
    INTEGER :: NODE,NPAR ! Number of ODEs, number of parameters.
    INTEGER :: LEFTBC ! Number of left boundary conditions.
    INTEGER :: NPTS ! Number of meshpoints. 
    INTEGER :: MXNSUB ! Maximum number of subintervals allowed.

    INTEGER :: NPTSIN ! Number of points in input mesh.
    DOUBLE PRECISION, DIMENSION(SOLIN%NODE) :: YANEW,YBNEW ! Solution values at
                                                           ! new endpoints.
    DOUBLE PRECISION, DIMENSION(SOLIN%NODE) :: VTMP1,VTMP2 ! Work arrays.
     
    INTEGER :: DEGREE,IER ! Degree of approximation for extending the solution,
                          ! Error flag for array allocation.
    INTEGER :: MODE ! Records type of extension to be performed.
    DOUBLE PRECISION :: A,B ! Problem interval endpoints.

!   Temporary storage of array info from SOLIN.
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: X,PARAMETERS
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: Y
!----------------------------------------------------------------------------------
    write(6,*) 'Into Extend_2'
!   Assign scalar info from SOLIN structure to local variables.
    NODE = SOLIN%NODE
    NPAR = SOLIN%NPAR
    LEFTBC = SOLIN%LEFTBC
    NPTS = SOLIN%NPTS
    MXNSUB = SOLIN%MXNSUB

!   Allocate temporary storage for info from SOLIN array fields and assign
!   these fields to local arrays.
    ALLOCATE(X(NPTS),Y(NODE,NPTS),STAT=IER)
    CALL CHECK_STAT(IER)
    X = SOLIN%X
    Y = SOLIN%Y
    IF (NPAR > 0) THEN
      ALLOCATE(PARAMETERS(NPAR),STAT=IER)
      CALL CHECK_STAT(IER)
      PARAMETERS = SOLIN%PARAMETERS
    END IF

!   Determine required type of extension and set MODE and NPTS accordingly.
!   Save a copy of the input NPTS value for later use.
    NPTSIN = NPTS
    A = SOLIN%X(1)
    B = SOLIN%X(NPTS)
    
    MODE = 0
    IF (ANEW < A) THEN
      NPTS = NPTS + 1
      MODE = 1
    END IF
    IF (B < BNEW) THEN
      NPTS = NPTS + 1      
      IF (MODE == 0) THEN
        MODE = 2
      ELSE
        MODE = 3
      END IF    
    END IF

!  Set DEGREE based on value of ORDER. Default is constant extension but
!  set to linear if ORDER is 1. Uses SOLIN.
    DEGREE = 0
    IF (PRESENT(ORDER)) THEN
      IF (ORDER == 1) DEGREE = 1
    END IF

    IF (DEGREE == 0) THEN
      YANEW = SOLIN%Y(:,1)
      YBNEW = SOLIN%Y(:,NPTSIN)
    ELSE
!     Linear extrapolation  
      IF (ANEW < A) THEN
        CALL BVP_EVAL(SOLIN,A,VTMP1,VTMP2)
        YANEW = VTMP1 + (ANEW - A)*VTMP2  
      END IF    
      IF (B < BNEW) THEN
        CALL BVP_EVAL(SOLIN,B,VTMP1,VTMP2)
        YBNEW = VTMP1 + (BNEW - B)*VTMP2
      END IF
    END IF
!   In either case, YANEW and YBNEW are now set appropriately and
!   SOLIN is no longer needed.

!   Relevant info from SOLIN is now unpacked. Note that IWORK and WORK 
!   arrays of SOLIN will no longer be valid after the mesh and solution 
!   info is extended. Deallocate all array fields of SOLIN.
    DEALLOCATE(SOLIN%X,SOLIN%Y,SOLIN%IWORK,SOLIN%WORK,STAT=IER)
    CALL CHECK_STAT(IER)
    IF (NPAR > 0) THEN
      DEALLOCATE(SOLIN%PARAMETERS,STAT=IER)
      CALL CHECK_STAT(IER)
    END IF


!   Pack solution info into SOLOUT structure; allocate array fields as needed.
    SOLOUT%INFO = 0
    SOLOUT%NODE = NODE
    SOLOUT%NPAR = NPAR
    SOLOUT%LEFTBC = LEFTBC
    SOLOUT%NPTS = NPTS

    IF (PRESENT(MAX_NUM_SUBINTERVALS)) THEN
      SOLOUT%MXNSUB = MAX_NUM_SUBINTERVALS
    ELSE
      SOLOUT%MXNSUB = MXNSUB
    ENDIF

    IF (NPAR > 0) THEN
      ALLOCATE(SOLOUT%PARAMETERS(NPAR),STAT=IER)
      CALL CHECK_STAT(IER)
      SOLOUT%PARAMETERS = PARAMETERS
    END IF
    
!   If P is present, assign parameter values to SOLOUT field. Verify that
!   the number of parameters is not changed.
    IF (PRESENT(P)) THEN
      IF (NPAR /= SIZE(P)) THEN
        PRINT *,'Cannot change the number of unknown parameters.'
        STOP
      END IF
      SOLOUT%PARAMETERS = P
    END IF

!   Allocate arrays for mesh and solution in SOLOUT.     
    ALLOCATE(SOLOUT%X(NPTS),SOLOUT%Y(NODE,NPTS),STAT=IER)
    CALL CHECK_STAT(IER)

!   Based on value of MODE, load mesh and solution fields appropriately.
    SELECT CASE(MODE)
      CASE (0) ! ANEW = A, B = BNEW
        SOLOUT%X = X
        SOLOUT%Y = Y
      CASE (1) ! ANEW < A, B = BNEW
        SOLOUT%X(1) = ANEW
        SOLOUT%X(2:NPTS) = X
        SOLOUT%Y(:,1) = YANEW
        SOLOUT%Y(:,2:NPTS) = Y
      CASE (2) ! ANEW = A, B < BNEW
        SOLOUT%X(1:NPTS-1) = X
        SOLOUT%X(NPTS) = BNEW
        SOLOUT%Y(:,1:NPTS-1) = Y
        SOLOUT%Y(:,NPTS) = YBNEW
      CASE (3) ! ANEW < A, B < BNEW
        SOLOUT%X(1) = ANEW
        SOLOUT%X(2:NPTS-1) = X
        SOLOUT%X(NPTS) = BNEW
        SOLOUT%Y(:,1) = YANEW
        SOLOUT%Y(:,2:NPTS-1) = Y
        SOLOUT%Y(:,NPTS) = YBNEW
    END SELECT
    
!   SOLOUT structure is now loaded. Deallocate temporary arrays.
    DEALLOCATE(X,Y,STAT=IER)
    CALL CHECK_STAT(IER)
    IF (NPAR > 0) THEN
      DEALLOCATE(PARAMETERS,STAT=IER)
      CALL CHECK_STAT(IER)
    END IF
    
  END FUNCTION EXTEND_2 

!=====================end of generic BVP_EXTEND==========================
  


  SUBROUTINE BVP_TERMINATE(SOL)
! BVP_TERMINATE deallocates the allocateable array fields of SOL, namely
! X, Y, PARAMETRERS, IWORK, and WORK.
!--------------------------------------------------------------------------
!     CALLED BY: User's program
!     CALLS: CHECK_STAT
!--------------------------------------------------------------------------
!   Input/output argument:
    TYPE(BVP_SOL) :: SOL

!   Local variables:    
    INTEGER :: IER
!----------------------------------------------------------------------------
    DEALLOCATE(SOL%X,SOL%Y,SOL%IWORK,SOL%WORK,STAT=IER)
    CALL CHECK_STAT(IER)

    IF (SOL%NPAR > 0) THEN
      DEALLOCATE(SOL%PARAMETERS,STAT=IER)
      CALL CHECK_STAT(IER) 
    END IF

  END SUBROUTINE BVP_TERMINATE

  END MODULE BVP_M
