! COLNEWSC - December 17, 2023 
!
!************************************************************************************************************************************************
!
! Authors: Dr.Paul Muir and Mark Adams, Saint Mary's University,Halifax,Canada
!
!
! This is a BETA version of COLNEWSC. Please report bugs to Mark Adams at adams.mark.alex@gmail.com with COLNEWSC
!in the subject line.
!
! COLNEWSC is a new software package that uses much of the COLNEW source code and, like COLNEW, is used to find numerical solutions for 
! boundary value ordinary differential equations.
!
! The main feature of the new software package is the use of a superconvergent interpolant as described in 
! Muir and Enright, Superconvergent interpolants for the collocation solution of boundary value ordinary differential equations,
! 1999. The superconvergent technology is only employed for K=2,3, or 4 collocation points. For K=1,5,6,7, the COLNEWSC executes identical to COLNEW.
! 
!
! Subroutines most affected are NEWMSH and ERRCHK. The source code which implements the superconvergent interpolant technology can
! be found at the end of this file. The NEWMSH routine performs mesh selection based on error "rough" error estimates computed
! in the NEWMSH routines otherwise it performs a simple mesh doubling. The error estimates in the ERRCHK routine are used to check for termination.
!
! Other new features in COLNEWSC include the capability to approximate the derivatives information and make DFSUB and DGSUB optional
! routines. DFSUB and DGSUB store the analytic derivative information.
!
! COLNEW is written in Fortran 77 whereas COLNEWSC has been adapted to FORTRAN 90 and uses some of FORTRAN 90 features. All of the
! following source code is contained within the COLNEW_M module.
!
! COLNEW represents the continuous approximate solution, i.e., the collocation solution, in terms of a piecewise polynomial basis
! of degree 'p' defined on a mesh that partitions the problem domain [a,b]. The unknown basis coefficients are obtained from the solution of a nonlinear
! system which consists of the boundary conditions together with a set of collocation equations. The collocation equations are
! obtained by requiring that the collocation solution satisfy the ODE system at the images of the set of 'p-1' Gauss points mapped onto each subinterval of the
! mesh. The nonlinear system is solved using a modified Newton iteration (based on quasilinearization) to obtain the basis
! coefficients. The computation is proportional to, 'N', the number of mesh subintervals.
!
! COLNEW uses two error estimation schemes - a primary, high quality error estimation scheme for checking the termination condition
! and a secondary "rough" error estimation scheme for mesh refinement. The primary error estimate is computed in the ERRCHK routine
! and  obtained using Richardson extrapolation. This means that the error estimate is obtained by computing the difference between the collocation solution 
! computed on a given mesh and a second collocation solution computed on a second mesh
! obtained by halving each subinterval of the first mesh. The error estimates  are obtained by computing the
! difference between the collocation solutions at several points within the subinterval. If the error estimates satisfy the user
! tolerance, the code terminates successfully. The ERRCHK routine can only be called after a mesh doubling.
!
! COLNEWSC, for 'k' collocation points =2,3, or 4, begins by computing primary error estimates only at the superconvergent  meshpoints. Once these meshpoint error estimates satisfy the user
! tolerance, then in the last part of the computation, COLNEWSC constructs two superconvergent interpolant solution approximations
! (SCSAs), one for the current mesh and one for the previous mesh sampled at several points within each subinterval. If this error estimate satisfies the user
! tolerance, the computation terminates and the SCSA is returned to the user as the continuous solution approximation. In most
! cases, this error estimate for the SCSA does in fact satisfy the user tolerance. If this is not the case, the current mesh is
! doubled and then a new collocation solution and corresponding SCSA are computed on the new mesh. This process is repeated until
! the user tolerance is satisfied. 
!
! The secondary error estimation scheme in the NEWMSH routine is used to predict the
! required number of mesh points of the next mesh and also determine if the locations of the points of the next mesh should be
! determined so that the error estimates are equidistributed over the subintervals of this next mesh. If this is not deemed to be
! worthwhile, the next mesh is obtained simply by doubling the current mesh. 
!
! The COLNEW_M module below is used by an module called COLNEW_WRAPPER which is stored in a separate file. First a driver (test) program, stored inside test.f90,
! uses the module COLNEW_WRAPPER inside the colnew_wrapper.f90 file and subsequently the module COLNEW_WRAPPER uses the COLNEW_M module inside this file. The main purpose of the colnew_wrapper.f90 file is to
! hide information from the user and simply the user interface in the driver program.
! 
!*********************************************************************************************************************************************************


MODULE COLNEW_M

! All the subroutines in this module are declared private except for 'colnew' and 'appsln'  below which are available to the colnew_wrapper_module.

PRIVATE 

PUBLIC:: COLNEW, APPSLN ! 

!
! Global variables
!-------------------
!
! 'switch' = 0 - implies primary error estimates are computed only at the superconvergent meshpoints. This is the first part of the computation.
! 
! 'switch' = 1 - implies primary error estimates are computed using two superconvergent interpolant solution approximations (SCSAs),  one for the current mesh and one 
!                for the previous mesh sampled at several points within each subinterval. This is the last part of the computation.
!

INTEGER :: SWITCH=0

! 'y_meshnew' and 'y_meshold' - solution approximations used to compute primary error estimates only at the meshpoints when 'switch' = 0.
!
! 'y_meshnew' - solution approximations for the current mesh.
!
! 'y_meshold' - solution approximations for the previous mesh.
!
! 'ykstore' and 'yknew' - integer counters to access parts of the 'y_meshnew' and 'y_meshold' arrays
!


DOUBLE PRECISION :: Y_MESHNEW(1000000),Y_MESHOLD(1000000)

INTEGER :: YKSTORE,YKNEW

! 'zold', 'dmzold', and 'meshold' - store the values for 'z','dmz', and 'xi' which are used in the primary error estimation when 'switch' = 1.
!
! 'zold' - solution approximations for the previous mesh.
!
! 'dmzold' - derivative approximations for the previous mesh.
!
! 'meshold' - meshpoint values for the previous mesh. 

DOUBLE PRECISION :: ZOLD(1000000),DMZOLD(1000000),MESHOLD(100000)

!
!

CONTAINS


!**********************************************************************
!  this package solves boundary value problems for
!  ordinary differential equations, as described below.
!
!  COLNEW is a modification of the package COLSYS by ascher,
!  christiansen and russell [1]. It incorporates a new basis
!  representation replacing b-splines, and improvements for
!  the linear and nonlinear algebraic equation solvers.
!  the package can be referenced as either COLNEW or COLSYS.
!**********************************************************************
!----------------------------------------------------------------------
!                            p a r t  1
!        main storage allocation and program control subroutines
!----------------------------------------------------------------------
!
      SUBROUTINE COLNEW (NCOMP, M, ALEFT, ARIGHT, ZETA, IPAR, LTOL, &
                         TOL, FIXPNT, ISPACE, FSPACE, IFLAG, &
                         FSUB, GSUB, GUESS, DFSUB, DGSUB)
!
!
!**********************************************************************
!
!     written by
!                  u. ascher,
!                            department of computer science,
!                            university of british columbia,
!                            vancouver, b. c., canada   v6t 1w5
!                  g. bader,
!                            institut f. angewandte mathematik
!                            university of heidelberg
!                            im neuenheimer feld 294
!                            d-6900 heidelberg 1
!
!**********************************************************************
!
!     purpose
!
!     this package solves a multi-point boundary value
!     problem for a mixed order system of ode-s given by
!
!          (m(i))
!         u       =  f  ( x; z(u(x)) )      i = 1, ... ,ncomp
!          i          i
!
!                                          aleft .lt. x .lt. aright,
!
!
!         g  ( zeta(j); z(u(zeta(j))) ) = 0   j = 1, ... ,mstar
!          j
!                                    mstar = m(1)+m(2)+...+m(ncomp),
!
!
!         where                          t
!               u = (u , u , ... ,u     )  is the exact solution vector
!                     1   2        ncomp
!
!                (mi)
!               u     is the mi=m(i) th  derivative of u
!                i                                      i
!
!                                  (1)        (m1-1)       (mncomp-1)
!               z(u(x)) = ( u (x),u  (x),...,u    (x),...,u      (x) )
!                            1     1          1            ncomp
!
!                f (x,z(u))   is a (generally) nonlinear function of
!                 i
!                             z(u)=z(u(x)).
!
!                g (zeta(j);z(u))  is a (generally) nonlinear function
!                 j
!                               used to represent a boundary condition.
!
!         the boundary points satisfy
!               aleft .le. zeta(1) .le. .. .le. zeta(mstar) .le. aright
!
!         the orders mi of the differential equations satisfy
!                            1 .le. m(i) .le. 4.
!
!
!**********************************************************************
!
!     method
!
!        the method used to approximate the solution u is
!     collocation at gaussian points, requiring m(i)-1 continuous
!     derivatives in the i-th component, i = 1, ..., ncomp.
!     here, k is the number of collocation points (stages) per
!     subinterval and is chosen such that k .ge. max m(i).
!     a runge-kutta-monomial solution representation is utilized.
!
!     references
!
!     [1] u. ascher, j. christiansen and r.d. russell,
!         collocation software for boundary-value odes,
!         acm trans. math software 7 (1981), 209-222.
!         this paper contains EXAMPLES where use of the code
!         is demonstrated.
!
!     [2] g. bader and u. ascher,
!         a new basis implementation for a mixed order
!         boundary value ode solver,
!         siam j. scient. stat. comput. (1987).
!
!     [3] u. ascher, j. christiansen and r.d. russell,
!         a collocation solver for mixed order
!         systems of boundary value problems,
!         math. comp. 33 (1979), 659-679.
!
!     [4] u. ascher, j. christiansen and r.d. russell,
!         colsys - a collocation code for boundary
!         value problems,
!         lecture notes comp.sc. 76, springer verlag,
!         b. childs et. al. (eds.) (1979), 164-185.
!
!     [5] c. deboor and r. weiss,
!         solveblok: a package for solving almost block diagonal
!         linear systems,
!         acm trans. math. software 6 (1980), 80-87.
!
!**********************************************************************
!
!     ***************     input to colnew     ***************
!
!     variables
!
!     ncomp - no. of differential equations   (ncomp .le. 20)
!
!     m(j) - order of the j-th differential equation
!            ( mstar = m(1) + ... + m(ncomp) .le. 40 )
!
!     aleft - left end of interval
!
!     aright - right end of interval
!
!     zeta(j) - j-th side condition point (boundary point). must
!               have  zeta(j) .le. zeta(j+1). all side condition
!               points must be mesh points in all meshes used,
!               see description of ipar(11) and fixpnt below.
!
!     ipar - an integer array dimensioned at least 11.
!            a list of the parameters in ipar and their meaning follows
!            some parameters are renamed in colnew; these new names are
!            given in parentheses.
!
!     ipar(1)     ( = nonlin )
!             = 0 if the problem is linear
!             = 1 if the problem is nonlinear
!
!     ipar(2) = no. of collocation points per subinterval  (= k )
!               where max m(i) .le.  k .le. 7 . if ipar(2)=0 then
!               colnew sets  k = max ( max m(i)+1, 5-max m(i) )
!
!     ipar(3) = no. of subintervals in the initial mesh  ( = n ).
!               if ipar(3) = 0 then colnew arbitrarily sets n = 5.
!
!     ipar(4) = no. of solution and derivative tolerances.  ( = ntol )
!               we require  0 .lt. ntol .le. mstar.
!
!     ipar(5) = dimension of fspace.     ( = ndimf )
!
!     ipar(6) = dimension of ispace.     ( = ndimi )
!
!     ipar(7) -  output control ( = iprint )
!              = -1 for full diagnostic printout
!              = 0 for selected printout
!              = 1 for no printout
!
!     ipar(8)     ( = iread )
!             = 0 causes colnew to generate a uniform initial mesh.
!             = 1 if the initial mesh is provided by the user.  it
!                 is defined in fspace as follows:  the mesh
!                 aleft=x(1).lt.x(2).lt. ... .lt.x(n).lt.x(n+1)=aright
!                 will occupy  fspace(1), ..., fspace(n+1). the
!                 user needs to supply only the interior mesh
!                 points  fspace(j) = x(j), j = 2, ..., n.
!             = 2 if the initial mesh is supplied by the user
!                 as with ipar(8)=1, and in addition no adaptive
!                 mesh selection is to be done.
!
!     ipar(9)     ( = iguess )
!             = 0 if no initial guess for the solution is
!                 provided.
!             = 1 if an initial guess is provided by the user
!                 in subroutine  guess.
!             = 2 if an initial mesh and approximate solution
!                 coefficients are provided by the user in  fspace.
!                 (the former and new mesh are the same).
!             = 3 if a former mesh and approximate solution
!                 coefficients are provided by the user in fspace,
!                 and the new mesh is to be taken twice as coarse;
!                 i.e.,every second point from the former mesh.
!             = 4 if in addition to a former initial mesh and
!                 approximate solution coefficients, a new mesh
!                 is provided in fspace as well.
!                 (see description of output for further details
!                 on iguess = 2, 3, and 4.)
!
!     ipar(10)= 0 if the problem is regular
!             = 1 if the first relax factor is =rstart, and the
!                 nonlinear iteration does not rely on past covergence
!                 (use for an extra sensitive nonlinear problem only).
!             = 2 if we are to return immediately upon  (a) two
!                 successive nonconvergences, or  (b) after obtaining
!                 error estimate for the first time.
!
!     ipar(11)= no. of fixed points in the mesh other than aleft
!               and aright. ( = nfxpnt , the dimension of fixpnt)
!               the code requires that all side condition points
!               other than aleft and aright (see description of
!               zeta ) be included as fixed points in fixpnt.
!
!     ltol  -  an array of dimension ipar(4). ltol(j) = l  specifies
!              that the j-th tolerance in  tol  controls the error
!              in the l-th component of z(u).   also require that
!              1.le.ltol(1).lt.ltol(2).lt. ... .lt.ltol(ntol).le.mstar
!
!     tol    - an array of dimension ipar(4). tol(j) is the
!              error tolerance on the ltol(j) -th component
!              of z(u). thus, the code attempts to satisfy
!              for j=1,...,ntol  on each subinterval
!              abs(z(v)-z(u))       .le. tol(j)*abs(z(u))       +tol(j)
!                            ltol(j)                     ltol(j)
!
!              if v(x) is the approximate solution vector.
!
!     fixpnt - an array of dimension ipar(11).   it contains
!              the points, other than aleft and aright, which
!              are to be included in every mesh.
!
!     ispace - an integer work array of dimension ipar(6).
!              its size provides a constraint on nmax,
!              the maximum number of subintervals. choose
!              ipar(6) according to the formula
!                      ipar(6)  .ge.  nmax*nsizei
!                where
!                      nsizei = 3 + kdm
!                with
!                      kdm = kd + mstar  ;  kd = k * ncomp ;
!                      nrec = no. of right end boundary conditions.
!
!
!     fspace - a real work array of dimension ipar(5).
!              its size provides a constraint on nmax.
!              choose ipar(5) according to the formula
!                      ipar(5)  .ge.  nmax*nsizef
!                where
!                      nsizef = 4 + 3 * mstar + (5+kd) * kdm +
!                              (2*mstar-nrec) * 2*mstar.
!
!
!     iflag - the mode of return from colnew.
!           = 1 for normal return
!           = 0 if the collocation matrix is singular.
!           =-1 if the expected no. of subintervals exceeds storage
!               specifications.
!           =-2 if the nonlinear iteration has not converged.
!           =-3 if there is an input data error.
!
!
!**********************************************************************
!
!     *************    user supplied subroutines   *************
!
!
!     the following subroutines must be declared external in the
!     main program which calls colnew.
!
!
!     fsub  - name of subroutine for evaluating f(x,z(u(x))) =
!                            t
!             (f ,...,f     )  at a point x in (aleft,aright).  it
!               1      ncomp
!             should have the heading
!
!                       subroutine fsub (x , z , f)
!
!             where f is the vector containing the value of fi(x,z(u))
!             in the i-th component and                            t
!                                       z(u(x))=(z(1),...,z(mstar))
!             is defined as above under  purpose .
!
!
!     dfsub - name of subroutine for evaluating the jacobian of
!             f(x,z(u)) at a point x.  it should have the heading
!
!                       subroutine dfsub (x , z , df)
!
!             where z(u(x)) is defined as for fsub and the (ncomp) by
!             (mstar) array df should be filled by the partial deriv-
!             atives of f, viz, for a particular call one calculates
!                                df(i,j) = dfi / dzj, i=1,...,ncomp
!                                                     j=1,...,mstar.
!
!
!     gsub  - name of subroutine for evaluating the i-th component of
!             g(x,z(u(x))) = g (zeta(i),z(u(zeta(i)))) at a point x =
!                             i
!             zeta(i) where 1.le.i.le.mstar. it should have the heading
!
!                       subroutine gsub (i , z , g)
!
!             where z(u) is as for fsub, and i and g=g  are as above.
!                                                     i
!             note that in contrast to f in  fsub , here
!             only one value per call is returned in g.
!
!
!     dgsub - name of subroutine for evaluating the i-th row of
!             the jacobian of g(x,u(x)).  it should have the heading
!
!                       subroutine dgsub (i , z , dg)
!
!             where z(u) is as for fsub, i as for gsub and the mstar-
!             vector dg should be filled with the partial derivatives
!             of g, viz, for a particular call one calculates
!                   dg(i,j) = dgi / dzj      j=1,...,mstar.
!
!
!     guess - name of subroutine to evaluate the initial
!             approximation for  z(u(x)) and for dmval(u(x))= vector
!             of the mj-th derivatives of u(x). it should have the
!             heading
!
!                       subroutine guess (x , z , dmval)
!
!             note that this subroutine is needed only if using
!             ipar(9) = 1, and then all  mstar  components of z
!             and  ncomp  components of  dmval  should be specified
!             for any x,  aleft .le. x .le. aright .
!
!
!**********************************************************************
!
!     ************   use of output from colnew   ************
!
!                 ***   solution evaluation   ***
!
!     on return from colnew, the arrays fspace and ispace
!     contain information specifying the approximate solution.
!     the user can produce the solution vector  z( u(x) )  at
!     any point x, aleft .le. x .le. aright, by the statement,
!
!           call appsln (x, z, fspace, ispace)
!
!     when saving the coefficients for later reference, only
!     ispace(1),...,ispace(7+ncomp)    and
!     fspace(1),...,fspace(ispace(7))    need to be saved as
!     these are the quantities used by appsln.
!
!
!                 ***   simple continuation   ***
!
!
!     a formerly obtained solution can easily be used as the
!     first approximation for the nonlinear iteration for a
!     new problem by setting   (iguess =) ipar(9) = 2, 3 or 4.
!
!     if the former solution has just been obtained then the
!     values needed to define the first approximation are
!     already in ispace and fspace.
!     alternatively, if the former solution was obtained in a
!     previous run and its coefficients were saved then those
!     coefficients must be put back into
!     ispace(1),..., ispace(7+ncomp)    and
!     fspace(1),..., fspace(ispace(7)).
!
!     for ipar(9) = 2 or 3 set ipar(3) = ispace(1) ( = the
!     size of the previous mesh ).
!
!     for ipar(9) = 4 the user specifies a new mesh of n subintervals
!     as follows.
!     the values in  fspace(1),...,fspace(ispace(7))  have to be
!     shifted by n+1 locations to  fspace(n+2),..,fspace(ispace(7)+n+1)
!     and the new mesh is then specified in fspace(1),..., fspace(n+1).
!     also set ipar(3) = n.
!
!
!**********************************************************************
!
!     ***************      package subroutines      ***************
!
!     the following description gives a brief overview of how the
!     procedure is broken down into the subroutines which make up
!     the package called  colnew . for further details the
!     user should refer to documentation in the various subroutines
!     and to the references cited above.
!
!     the subroutines fall into four groups:
!
! part 1 - the main storage allocation and program control subr
!
!     colnew - tests input values, does initialization and breaks up
!              the work areas, fspace and ispace, into the arrays
!              used by the program.
!     colsys - another name for colnew
!
!     contrl - is the actual driver of the package. this routine
!              contains the strategy for nonlinear equation solving.
!
!     skale  - provides scaling for the control
!              of convergence in the nonlinear iteration.
!
!
! part 2 - mesh selection and error estimation subroutines
!
!     consts - is called once by  colnew  to initialize constants
!              which are used for error estimation and mesh selection.
!
!     newmsh - generates meshes. it contains the test to decide
!              whether or not to redistribute a mesh.
!
!     errchk - produces error estimates and checks against the
!              tolerances at each subinterval
!
!
! part 3 - collocation system set-up subroutines
!
!     lsyslv - controls the set-up and solution of the linear
!              algebraic systems of collocation equations which
!              arise at each newton iteration.
!
!              Other subroutines used by lsyslv -
!              
!              1) dmzsol - computes dmz in a blockwise manner
!              2) fcblok - (which calls factrb and shiftb)
!
!                 fcblok  supervises the plu factorization with pivoting of
!                 scaled rows of the almost block diagonal matrix stored in the
!                 arrays  bloks  and  integs .
!
!                 factrb = subprogram which carries out steps 1,...,last of gauss
!                 elimination (with pivoting) for an individual block.
!
!                 shiftb = subprogram which shifts the remaining rows to the top of
!                 the next block
!
!              3) sbblok (which calls subfor and subbak)
!
!                 supervises the solution (by forward and backward substitution) of
!                 the linear system  a*x = b  for x, with the plu factorization of
!                 a  already generated in  fcblok .  individual blocks of
!                 equations are solved via  subfor  and  subbak .
!
!     gderiv - is used by lsyslv to set up the equation associated
!              with a side condition point.
!
!     vwblok - is used by lsyslv to set up the equation(s) associated
!              with a collocation point.
!
!     gblock - is used by lsyslv to construct a block of the global
!              collocation matrix or the corresponding right hand
!              side.
!
!
! part 4 - service subroutines
!
!     appsln - sets up a standard call to  approx .
!
!     approx - evaluates a piecewise polynomial solution.
!
!     rkbas  - evaluates the mesh independent runge-kutta basis
!
!     vmonde - solves a vandermonde system for given right hand
!              side
!
!     horder - evaluates the highest order derivatives of the
!              current collocation solution used for mesh refinement.
!
!
! part 5 - linear algebra  subroutines
!
!     to solve the global linear systems of collocation equations
!     constructed in part 3,  colnew  uses a column oriented version
!     of the package  solveblok originally due to de boor and weiss.
!
!     to solve the linear systems for static parameter condensation
!     in each block of the collocation equations, the linpack
!     routines  dgefa and  dgesl  are included. but these
!     may be replaced when solving problems on vector processors
!     or when solving large scale sparse jacobian problems.
!
!----------------------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION M(1), ZETA(1), IPAR(11), LTOL(1), TOL(1), DUMMY(1), &
                FIXPNT(1),ISPACE(*), FSPACE(*)


      COMMON /COLOUT/ PRECIS, IOUT, IPRINT
      COMMON /COLLOC/ RHO(7), COEF(49)
      COMMON /COLORD/ K, NC, MSTAR, KD, MMAX, MT(20)
      COMMON /COLAPR/ N, NOLD, NMAX, NZ, NDMZ
      COMMON /COLMSH/ MSHFLG, MSHNUM, MSHLMT, MSHALT
      COMMON /COLSID/ TZETA(40), TLEFT, TRIGHT, IZETA, IDUM
      COMMON /COLNLN/ NONLIN, ITER, LIMIT, ICARE, IGUESS
      COMMON /COLEST/ TTL(40), WGTMSH(40), WGTERR(40), TOLIN(40), &
                      ROOT(40), JTOL(40), LTTOL(40), NTOL
!
      EXTERNAL FSUB, DFSUB, GSUB, DGSUB, GUESS

!     The 'dfsub' and 'dgsub' subroutines provided by the user are optional
!

      OPTIONAL :: DFSUB, DGSUB
!
!*********************************************************************
!
!     the actual subroutine colnew serves as an interface with
!     the package of subroutines referred to collectively as
!     colnew. the subroutine serves to test some of the input
!     parameters, rename some of the parameters (to make under-
!     standing of the coding easier), to do some initialization,
!     and to break the work areas fspace and ispace up into the
!     arrays needed by the program.
!
!**********************************************************************
!
!
!     Description of variables in all common blocks.
!     There are 9 common blocks in total and the only one this subroutine (colnew)
!     does not use is colbas
!
!     colout- 1) Precis is a constant = 100*macine unit roundoff
!             2) Iout is a constant = 6 for writing output to the screen
!             3) Iprint is equal to ipar(7) and therefore iprint
!         
!                = -1 for full diagnostic printout
!                = 0 for selected printout
!                = 1 for no printout
!      
!     colloc - 1) Rho are the k collocation points on (0,1)
!              2) Coef are the precomputed derivatives of the basis
!
!     colord - 1) K is the number of collocation points
!              2) Nc is equal to ncomp is the number of differential equations
!              3) Mstar is the sum of the orders of the differential equations
!                 which is equal to the number of solution and derivative tolerances
!                 which is also equal to the number of differential equations in an
!                 equivalent first order system.
!              4) Kd is the number of collocation points multiplied by the number of
!                 differential equations.
!              5) Mmax is the maxmium order of the set of differential equations
!              6) Mt is an array which stores all of the orders of each differential equation
!
!    colapr - 1) N is the number of subintervals in the current mesh
!             2) Nold is the number of subintervals in the previous mesh
!             3) Nmax is the maximum number of subintervals allowed
!             4) Nz is the number of components in the Z solution vector equal to
!                Mstar*(N+1)
!             5) Ndmz is equal to the number of collocation points multiplied by the number of differential
!                equations multiplied by the number of subintervals.
!                This is equal to the number of components of Dmz equal to Kd*N
!
!   colmsh - 1) mshflg = 1 implies the mesh is a halving of its former mesh (so an error estimate has been calculated)
!                      = 0 otherwise
!            2) mshnum is the number of mesh selections that have been performed for the given N.
!            3) mshlmt is the maximum number of mesh selections which are permitted for a given N before mesh halving.
!            4) mshalt is the nummber of consecutive times (plus 1) the mesh slection has alternatively halved and doubled n.
!               If mshalt is greater then or equal to mshlmt then contrl required that the current mesh be halved. 
!   colsid - 1) tzeta is an array equal to zeta(i) so
!
!               tzeta(j) - j-th side condition point (boundary point). must
!                          have  tzeta(j) .le. tzeta(j+1). all side condition
!                          points must be mesh points in all meshes used,
!                          see description of ipar(11) and fixpnt below. 
!
!               (this variable name is not used in computation)
!
!             2) tleft is equal to aleft the left boundary point
!                (this variable name is not used in computation)
!             3) tright is equal to aright the right boundary point
!                (this variable name is not used in computation)
!             4) izeta points to zeta (allows incrementation by 1)
!             5) idum equals izave equals izeta
!
!  colnln    1) nonlin is equal to ipar(1)
!                = 0 if the problem is linear
!                = 1 if the problem is nonlinear
!
!            2) iter is the number of iterations (counted only when Jacobian inversions are performed)
!            3) limit is the maximum number of newton iterations per mesh
!            4) icare = ipar(10) initially
!
!               = 0 if the problem is regular
!               = 1 if the first relax factor is =rstart, and the
!                 nonlinear iteration does not rely on past covergence
!                 (use for an extra sensitive nonlinear problem only).
!               = 2 if we are to return immediately upon  (a) two
!                 successive nonconvergences, or  (b) after obtaining
!                 error estimate for the first time.
!
!               Later in the program the meaning changes slightly to
!
!               = -1 no convergence occured for a reular problem
!               = 0 a regular problem
!               = 1 a sensitive problem
!               = 2 used for continuation ( see description of ipar(10)
!                 in colnew).
!             5) iguess = ipar(9) 
!                 = 0 if no initial guess for the solution is
!                 provided.
!                 = 1 if an initial guess is provided by the user
!                 in subroutine  guess.
!                 = 2 if an initial mesh and approximate solution
!                 coefficients are provided by the user in  fspace.
!                 (the former and new mesh are the same).
!                 = 3 if a former mesh and approximate solution
!                 coefficients are provided by the user in fspace,
!                 and the new mesh is to be taken twice as coarse;
!                 i.e.,every second point from the former mesh.
!                 = 4 if in addition to a former initial mesh and
!                 approximate solution coefficients, a new mesh
!                 is provided in fspace as well.
!                 (see description of output for further details
!                 on iguess = 2, 3, and 4. see also simple continuation documentation above)
!
!     colest - 1) ttl is equal to the error tolerances 
!              (variable name not actually used in computation)
!              2)wgtmsh is gotten from the values of cnsts2 which in turn are
!              the constants in the theoretical expression for the
!              errors. the quantities in wgtmsh are 10x the values
!              in cnsts2 so that the mesh selection algorithm
!              is aiming for errors .1x as large as the user
!              requested tolerances.i
!              (Note:
!    
!                  arrays assigned during compilation:
!                  cnsts1 - weights for extrapolation error estimate
!                  cnsts2 - weights for mesh selection
!                  (the above weights come from the theoretical form for
!                   the collocation error -- see [3])
!
!             3) wgterr - the particular values of cnsts1 used for current run
!              (depending on k, m)
!           
!             4) tolin equals the error tolerances
!             5) root   - reciprocals of expected rates of convergence of compo-
!                         nents of z(j) for which tolerances are specified
!             6) jtol   - components of differential system to which tolerances
!                refer (viz, if ltol(i) refers to a derivative of u(j),
!                then jtol(i)=j)
!    
!             7) lttol = ltol (i) 
!
!                 ltol  -  an array of dimension ipar(4). ltol(j) = l  specifies
!                 that the j-th tolerance in  tol  controls the error
!                 in the l-th component of z(u).   also require that
!                  1.le.ltol(1).lt.ltol(2).lt. ... .lt.ltol(ntol).le.mstar
!             8) ntol is the number of solution and derivative tolerances
!
!       colbas - 1) B=
!                   = rkb - the runge kutta basis (0th to (m-1)th derivatives)
!                   computed in rkbas at S (the relative position for which the basis is to be evaluated)
!                   (0.le.S.le.1)
!                     
!                2) acol is the runge kutta coefficients at the collocation points (rho) = S
!                3) asave is the runge-kutta coefficients at (1/6,1/3,2/3,5/6) =S used in
!                   newmesh and errchk 
!   
!----------------------------------------------------------------------------------------------------------------------
!
!     Called by subroutines : colnew_solver (inside colnew_wrapper module) 
!               
!     Calls to subroutines : consts, newmsh, contrl
!
!-----------------------------------------------------------------------------------------------------------------------
!
!...  specify machine dependent output unit  iout  and compute machine
!...  dependent constant  precis = 100 * machine unit roundoff
!
!
!     initialize 'switch' to zero
!
      SWITCH = 0
!  
!
      IF ( IPAR(7) .LE. 0 )  WRITE(6,99)
  99  FORMAT(//,33H VERSION *COLNEW* OF COLSYS .    ,//)
!
!     Calculate the value for 100 * machine unit roundoff
!
      IOUT = 6
      PRECIS = 1.D0
   10 PRECIS = PRECIS / 2.D0
      PRECP1 = PRECIS + 1.D0
      IF ( PRECP1 .GT. 1.D0 )                       GO TO 10
      PRECIS = PRECIS * 100.D0
!
!...  in case incorrect input data is detected, the program returns
!...  immediately with iflag=-3.
!
!     Only currently enough memory for 20 equations
!
      IFLAG = -3
      IF ( NCOMP .LT. 1 .OR. NCOMP .GT. 20 )        RETURN

!     The order of the differntial equations must be 4 or less.

      DO 20 I=1,NCOMP
         IF ( M(I) .LT. 1 .OR. M(I) .GT. 4 )        RETURN
   20 CONTINUE
!
!...  rename some of the parameters and set default values.
!
      NONLIN = IPAR(1)
      K = IPAR(2)
      N = IPAR(3)
      IF ( N .EQ. 0 )  N = 5


!     Iread is a local variable equal to ipar(8) such that:
!             = 0 causes colnew to generate a uniform initial mesh.
!             = 1 if the initial mesh is provided by the user.  it
!                 is defined in fspace as follows:  the mesh
!                 aleft=x(1).lt.x(2).lt. ... .lt.x(n).lt.x(n+1)=aright
!                 will occupy  fspace(1), ..., fspace(n+1). the
!                 user needs to supply only the interior mesh
!                 points  fspace(j) = x(j), j = 2, ..., n.
!             = 2 if the initial mesh is supplied by the user
!                 as with ipar(8)=1, and in addition no adaptive
!                 mesh selection is to be done.


      IREAD = IPAR(8)
      IGUESS = IPAR(9)

!     If the equations are linear and then no initial guess is required even if it supplied by the user subroutine.

      IF ( NONLIN .EQ. 0 .AND. IGUESS .EQ. 1 )  IGUESS = 0

!     If the new guess is obtained from the former mesh (simple continuation) then iread =1 implies the corresponding mesh
!     has to be supplied as well from fspace. Iread =0 does not make sense because we do not want a uniform initial mesh automatically generated.

      IF ( IGUESS .GE. 2 .AND. IREAD .EQ. 0 )   IREAD = 1
      ICARE = IPAR(10)
      NTOL = IPAR(4)

!     Ndimf and Ndimi are local variables which store the dimensions of fspace and ispace respectively.

      NDIMF = IPAR(5)
      NDIMI = IPAR(6)
!
!     Nfixpnt is a local variable equal to the number of extra fixpoints.

      NFXPNT = IPAR(11)
      IPRINT = IPAR(7)
      MSTAR = 0
      MMAX = 0
      DO  30 I = 1, NCOMP
         MMAX = MAX0 ( MMAX, M(I) )
         MSTAR = MSTAR + M(I)
         MT(I) = M(I)
   30 CONTINUE
      IF ( K .EQ. 0 )   K = MAX0( MMAX + 1 , 5 - MMAX )
      DO 40 I = 1, MSTAR
   40 TZETA(I) = ZETA(I)
      DO 50 I = 1, NTOL
         LTTOL(I) = LTOL(I)
   50 TOLIN(I) = TOL(I)
      TLEFT = ALEFT
      TRIGHT = ARIGHT
      NC = NCOMP
      KD = K * NCOMP
!
!...  Print the input data for checking.
!
      IF ( IPRINT .GT. -1 )                         GO TO 80
      IF ( NONLIN .GT. 0 )                          GO TO 60
      WRITE (IOUT,260) NCOMP, (M(IP), IP=1,NCOMP)
      GO TO 70
   60 WRITE (IOUT,270) NCOMP, (M(IP), IP=1,NCOMP)
   70 WRITE (IOUT,280) (ZETA(IP), IP=1,MSTAR)
      IF ( NFXPNT .GT. 0 ) &
         WRITE (IOUT,340) NFXPNT, (FIXPNT(IP), IP=1,NFXPNT)
      WRITE (IOUT,290) K
      WRITE (IOUT,300) (LTOL(IP), IP=1,NTOL)
      WRITE (IOUT,310) (TOL(IP), IP=1,NTOL)
      IF (IGUESS .GE. 2) WRITE (IOUT,320)
      IF (IREAD .EQ. 2) WRITE (IOUT,330)
   80 CONTINUE


!
!...  Check for correctness of data
!
      IF ( K .LT. 0 .OR. K .GT. 7 )                 RETURN
      IF ( N .LT. 0 )                               RETURN
      IF ( IREAD .LT. 0 .OR. IREAD .GT. 2 )         RETURN
      IF ( IGUESS .LT. 0 .OR. IGUESS .GT. 4 )       RETURN
      IF ( ICARE .LT. 0 .OR. ICARE .GT. 2 )         RETURN
      IF ( NTOL .LT. 0 .OR. NTOL .GT. MSTAR )       RETURN
      IF ( NFXPNT .LT. 0 )                          RETURN
      IF ( IPRINT .LT. (-1) .OR. IPRINT .GT. 1 )    RETURN
      IF ( MSTAR .LT. 0 .OR. MSTAR .GT. 40 )        RETURN
      IP = 1
      DO 100 I = 1, MSTAR
      IF ( DABS(ZETA(I) - ALEFT) .LT. PRECIS .OR. &
           DABS(ZETA(I) - ARIGHT) .LT. PRECIS )     GO TO 100
   90 IF ( IP .GT. NFXPNT )                         RETURN
        IF ( ZETA(I) - PRECIS .LT. FIXPNT(IP) )     GO TO 95
        IP = IP + 1
      GO TO 90
   95 IF ( ZETA(I) + PRECIS .LT. FIXPNT(IP) )       RETURN
  100 CONTINUE


!
!...  Set limits on iterations and initialize counters.
!...  Limit is a local variable = maximum number of newton iterations per mesh.
!
!
!     Only 3 mesh selections are permitted for a given N before mesh halving is done.

      MSHLMT = 3
!
!     The mesh is currently NOT a mesh halving of its former mesh (so an error estimate has NOT been calculated)

      MSHFLG = 0
!  
!     There has been only 1 inputted mesh selection done for the given N.

      MSHNUM = 1

!     There has been 1 consecutive time the selection process has halved the current mesh.

      MSHALT = 1


!...  Limit is a local variable = maximum number of newton iterations per mesh.

      LIMIT = 40
!
!...  Compute the maxium possible n for the given sizes of
!...  ispace  and  fspacei stored in the local variable nmax.
!
      NREC = 0
      DO 110 I = 1, MSTAR
           IB = MSTAR + 1 - I
           IF ( ZETA(IB) .GE. ARIGHT )  NREC = I
  110 CONTINUE
      NFIXI = MSTAR
      NSIZEI = 3 + KD + MSTAR
      NFIXF = NREC * (2*MSTAR) + 5 * MSTAR + 3
      NSIZEF = 4 + 3 * MSTAR + (KD+5) * (KD+MSTAR) + &
      (2*MSTAR-NREC) * 2*MSTAR

      NMAXF = (NDIMF - NFIXF) / NSIZEF
      NMAXI = (NDIMI - NFIXI) / NSIZEI


      IF ( IPRINT .LT. 1 )  WRITE(IOUT,350) NMAXF, NMAXI
      NMAX = MIN0( NMAXF, NMAXI )


      IF ( NMAX .LT. N )                            RETURN
      IF ( NMAX .LT. NFXPNT+1 )                     RETURN
      IF (NMAX .LT. 2*NFXPNT+2 .AND. IPRINT .LT. 1)  WRITE(IOUT,360)
!
!...  Generate pointers to break up  fspace  and  ispace .
!
!     LXI points to XI the current mesh of maxmium size NMAX + 1
!
!     LG points to G of maximum size 2*mstar*(Nmax*(2*Mstar-nrec)+nrec used in LSYSLV where, 
!
!         this routine controls the set up and solution of a linear
!      system of collocation equations.
!         the matrix  G  is cast into an almost block diagonal
!      form by an appropriate ordering of the columns and solved
!      using the package of de boor-weiss [5]. the matrix is composed
!      of n blocks. the i-th block has the size
!                  integs(1,i) * integs(2,i).
!      it contains in its last rows the linearized collocation
!      equations, condensed as described in [2],
!      and the linearized side conditions corresponding to
!      the i-th subinterval.  integs(3,i)  steps of gaussian
!      elimination are applied to it to achieve a  partial plu
!      decomposition.  the right hand side vector is put into  rhs
!      and the solution vector is returned in  delz and deldmz.
!
!     LXIOLD points XIOLD to the previous mesh of maximum size NMAX+1
!
!     LW and LV point to the blocks w and v of the collocation matrix corresponding to
!     an interior collocation point (before parameter condensation)
!
!     W is of maximum size KD**2*NMAX
!     V is of maximum size MSTAR*(KD)*NMAX
!  
!     LZ points to Z is the solution vector of maximum size Mstar*(Nmax+1)
!
!     LDMZ points to DMZ of maximum size KD*NMAX is the array of mj-th derivatives of the current
!     solution
!
!     LDELZ points to DELZ or delta z or the change in Z , part of the solution to the linear system
!     of maxmium size (MSTAR*(NMAX+1))
!
!     LDELDZ point to DELDMZ or delta dmz or change in dmz, part of the solution to the linear system of maxmium
!     size KD*NMAX
!
!     LDQZ is of maximum size (MSTAR*(NMAX+1)) 
!
!     LDQDMZ is of maximum size KD*NMAX
!
!     LRHS points to the right hand side of the linear system of maxmium size KD*NMAX+MSTAR
!
!     LVALST points to VALSTR of maxmium size 4*MSTAR*NMAX is assigned values need in errchk for the
!     error estimate
!
!     LSLOPE points to SLOPE is of maximum size NMAX is described in newmsh as

!            slope  - an approximate quantity to be equidistributed for
!                     mesh selection (see [3]), viz,
!                             .                        (k+mj)
!                     slope(i)=     max   (weight(l) *u      (xi(i)))
!                               1.le.l.le.ntol         j
!
!                     where j=jtol(l)
!
!     LACCUM points to ACCUM is of maximum size NMAX+1 where in newmsh it is described as 
!     accum(i) is the integral of the slope from aleft to xiold(i)
!
!     LSCL points to SCALE of maximum size MSTAR*(NMAX+1) is the scaling vector for Z
!
!     LDSCL points to DSKALE is the scaling vector for dmz
!
!     Now break up ispace which breaks up pointers to ipvtg, integs,ipvtw having to do with linear system solve
!     Ipvtg is of maximum size Mstar*(Nmax+1) and Ipvtw is of maximum size KD*Nmax


      LXI = 1
      LG = LXI + NMAX + 1
      LXIOLD = LG + 2*MSTAR * (NMAX * (2*MSTAR-NREC) + NREC)
      LW     = LXIOLD + NMAX + 1
      LV     = LW + KD**2 * NMAX
      LZ     = LV + MSTAR * KD * NMAX
      LDMZ   = LZ + MSTAR * (NMAX + 1)
      LDELZ  = LDMZ + KD * NMAX
      LDELDZ = LDELZ + MSTAR * (NMAX + 1)
      LDQZ   = LDELDZ + KD * NMAX
      LDQDMZ = LDQZ + MSTAR * (NMAX + 1)
      LRHS   = LDQDMZ + KD * NMAX
      LVALST = LRHS   + KD * NMAX + MSTAR
      LSLOPE = LVALST + 4 * MSTAR * NMAX
      LACCUM = LSLOPE + NMAX
      LSCL   = LACCUM + NMAX + 1
      LDSCL  = LSCL + MSTAR * (NMAX + 1)
      LPVTG = 1
      LPVTW = LPVTG + MSTAR * (NMAX + 1)
      LINTEG = LPVTW + KD * NMAX
!
!...  If  iguess .ge. 2, move  xiold, z, and  dmz  to their proper
!...  locations in  fspace. Save in fspace for the simple continuation process.
!
! Recall:


!                 ***   simple continuation   ***
!
!
!     a formerly obtained solution can easily be used as the
!     first approximation for the nonlinear iteration for a
!     new problem by setting   (iguess =) ipar(9) = 2, 3 or 4.
!
!     if the former solution has just been obtained then the
!     values needed to define the first approximation are
!     already in ispace and fspace.
!     alternatively, if the former solution was obtained in a
!     previous run and its coefficients were saved then those
!     coefficients must be put back into
!     ispace(1),..., ispace(7+ncomp)    and
!     fspace(1),..., fspace(ispace(7)).
!
!     for ipar(9) = 2 or 3 set ipar(3) = number of subintervals in the initial mesh = ispace(1) ( = the
!     size of the previous mesh ).
!
!     for ipar(9) = 4 the user specifies a new mesh of n subintervals
!     as follows.
!     the values in  fspace(1),...,fspace(ispace(7))  have to be
!     shifted by n+1 locations to  fspace(n+2),..,fspace(ispace(7)+n+1)
!     and the new mesh is then specified in fspace(1),..., fspace(n+1).
!     also set ipar(3) = n.
!
!
!**********************************************************************


!
      IF ( IGUESS .LT. 2 )                          GO TO 160
!     Go to 160 because no continuation is taking place.

      
      NOLD = N
      IF (IGUESS .EQ. 4)  NOLD = ISPACE(1) ! For ipar(9)=4 the user supplies the old value for N in ispace(1)
      ! NZ will now store the old value of the Z solution vector

      NZ = MSTAR * (NOLD + 1)

      ! NDMZ will now store the old value of dmz.

      NDMZ = KD * NOLD
       
      ! NP1 is a local variable storing the size of the new mesh.

      NP1 = N + 1
      IF ( IGUESS .EQ. 4 )  NP1 = NP1 + NOLD + 1 !For ipar(9)=4 store the number of mesh points in the previous mesh and the new mesh.
 

      DO 120 I=1,NZ
  120 FSPACE( LZ+I-1 )  =  FSPACE( NP1+I )


      IDMZ = NP1 + NZ
      DO 125 I=1,NDMZ
  125 FSPACE( LDMZ+I-1 )  =  FSPACE( IDMZ+I )
      NP1 = NOLD + 1
      IF ( IGUESS .EQ. 4 )                          GO TO 140
      DO 130 I=1,NP1
  130 FSPACE( LXIOLD+I-1 )  =  FSPACE( LXI+I-1 )
      GO TO 160
  140 DO 150 I=1,NP1
  150 FSPACE( LXIOLD+I-1 )  =  FSPACE( N+1+I )
  160 CONTINUE
!
!...  initialize collocation points, constants, mesh.
!
!     We have K so RHO and COEF are dependent on K.
!     where  Rho are the k collocation points on (0,1)
!     and Coef are the precomputed derivatives of the basis

      CALL CONSTS ( K, RHO, COEF )


!     Call newmsh with mode (the first argument) = 3 + Iread (see iread above = ipar(8) =0,1,2)
!     Most likely iread = 0 for first call to newmsh which mean to create a uniform initial mesh.
!     So 3 + Iread = 3+0 = 3.
!
!     Mode = 5,4,3 - deal mainly with definition of an initial
!                    mesh for the current boundary value problem

      CALL NEWMSH (3+IREAD, FSPACE(LXI), FSPACE(LXIOLD), DUMMY, &
                   DUMMY, DUMMY, DUMMY, DUMMY, NFXPNT, FIXPNT,FSUB)
!
!...  Determine first approximation, if the problem is nonlinear (most likely nonlinear).
!
!     Recall Iguess greater than 1 implies simple continuation and first approximation is provided in fspace.
!
      IF (IGUESS .GE. 2)                            GO TO 230
      NP1 = N + 1
      DO 210 I = 1, NP1

!     For input to control the old mesh has to be the same as the new mesh because if Iguess is greater than 2 the old or former  mesh is
!     is already provided.


  210 FSPACE( I + LXIOLD - 1 ) = FSPACE( I + LXI - 1 )
      NOLD = N

!     Recall iguess = 1 implies first approximation is provided in subroutine guess so then GOTO just before subroutine CONTRL next

      IF ( NONLIN .EQ. 0  .OR. IGUESS .EQ. 1 )      GO TO 230
!
!...  System provides first approximation of the solution if subroutine guess is not given.
!...  Simply choose z(j) = 0  for j=1,...,mstar.
!
      DO 220 I=1, NZ
  220 FSPACE( LZ-1+I ) = 0.D0
      DO 225 I=1, NDMZ
  225 FSPACE( LDMZ-1+I ) = 0.D0
  230 CONTINUE

!     Reset to no approximation given if using continuation for common block colnln. 

      IF (IGUESS .GE. 2)  IGUESS = 0

!     Call subroutine Contrl, the actual driver,  The nonlinear iteration
!     strategy is controlled here ( see [4] ). Upon convergence, errchk
!     is called to test for satisfaction of the requested tolerances.

      CALL CONTRL (FSPACE(LXI),FSPACE(LXIOLD),FSPACE(LZ),FSPACE(LDMZ), &
           FSPACE(LRHS),FSPACE(LDELZ),FSPACE(LDELDZ),FSPACE(LDQZ), &
           FSPACE(LDQDMZ),FSPACE(LG),FSPACE(LW),FSPACE(LV), &
           FSPACE(LVALST),FSPACE(LSLOPE),FSPACE(LSCL),FSPACE(LDSCL), &
           FSPACE(LACCUM),ISPACE(LPVTG),ISPACE(LINTEG),ISPACE(LPVTW), &
           NFXPNT,FIXPNT,IFLAG,FSUB,DFSUB,GSUB,DGSUB,GUESS )
!
!...  Prepare output
!
!     The size of the last mesh is returned.

      ISPACE(1) = N

!     The number of collocation points is returned.

      ISPACE(2) = K

!     The number of  of differential equations is returned

      ISPACE(3) = NCOMP

!     The number of components on the solution vector z for each mesh point is returned

      ISPACE(4) = MSTAR
!    
!     The maximum order in the set of differential equations is returned.

      ISPACE(5) = MMAX

!     Ispace(6) is used to locate the coef (see below) for later calls with user routine appsln.

      ISPACE(6) = NZ + NDMZ + N + 2
      K2 = K * K

!     Currently do not see where ispace(7) is used.

      ISPACE(7) = ISPACE(6) + K2 - 1

!     The orders of each differential equation are returned.

      DO 240 I = 1, NCOMP
  240 ISPACE(7+I) = M(I)

!     NZ is the number of components of the solution vector.
!     Return the values of the solution vector.

      DO 250 I = 1, NZ
  250 FSPACE( N+1+I ) = FSPACE( LZ-1+I )


!     NDMZ is the number of components in the dmz solution vector.
!     Return the dmz solution vector = to the mj-th derivatives of the solution vector

      IDMZ = N + 1 + NZ
      DO 255 I = 1, NDMZ
  255 FSPACE( IDMZ+I ) = FSPACE( LDMZ-1+I )
      IC = IDMZ + NDMZ

!     Return coef, the derivatives of the basis

      DO 258 I = 1, K2
  258 FSPACE( IC+I ) = COEF(I)
      RETURN
!----------------------------------------------------------------------
  260 FORMAT(/// 37H THE NUMBER OF (LINEAR) DIFF EQNS IS , I3/ 1X, &
             16HTHEIR ORDERS ARE, 20I3)
  270 FORMAT(/// 40H THE NUMBER OF (NONLINEAR) DIFF EQNS IS , I3/ 1X, &
             16HTHEIR ORDERS ARE, 20I3)
  280 FORMAT(27H SIDE CONDITION POINTS ZETA, 8F10.6, 4( / 27X, 8F10.6))
  290 FORMAT(37H NUMBER OF COLLOC PTS PER INTERVAL IS, I3)
  300 FORMAT(39H COMPONENTS OF Z REQUIRING TOLERANCES -,8(7X,I2,1X), &
             4(/38X,8I10))
  310 FORMAT(33H CORRESPONDING ERROR TOLERANCES -,6X,8D10.2, &
             4(/39X,8D10.2))
  320 FORMAT(44H INITIAL MESH(ES) AND Z,DMZ PROVIDED BY USER)
  330 FORMAT(27H NO ADAPTIVE MESH SELECTION)
  340 FORMAT(10H THERE ARE ,I5,27H FIXED POINTS IN THE MESH - , &
             10(6F10.6/))
  350 FORMAT(44H THE MAXIMUM NUMBER OF SUBINTERVALS IS MIN (, I4, &
             23H (ALLOWED FROM FSPACE),,I4, 24H (ALLOWED FROM ISPACE) ))
  360 FORMAT(/53H INSUFFICIENT SPACE TO DOUBLE MESH FOR ERROR ESTIMATE)


      END SUBROUTINE COLNEW



      SUBROUTINE CONTRL (XI, XIOLD, Z, DMZ, RHS, DELZ, DELDMZ,&
                 DQZ, DQDMZ, G, W, V, VALSTR, SLOPE, SCALE, DSCALE, &
                 ACCUM, IPVTG, INTEGS, IPVTW, NFXPNT, FIXPNT, IFLAG,&
                 FSUB, DFSUB, GSUB, DGSUB, GUESS )
!
!**********************************************************************
!
!   purpose
!     this subroutine is the actual driver.  the nonlinear iteration
!     strategy is controlled here ( see [4] ). upon convergence, errchk
!     is called to test for satisfaction of the requested tolerances.
!
!
!   variables
!
!     check  - maximum tolerance value, used as part of criteria for
!              checking for nonlinear iteration convergence
!     relax  - the relaxation factor for damped newton iteration
!     relmin - minimum allowable value for relax  (otherwise the
!              jacobian is considered singular).
!     rlxold - previous relax
!     rstart - initial value for relax when problem is sensitive
!     ifrz   - number of fixed jacobian iterations
!     lmtfrz - maximum value for ifrz before performing a reinversion
!     iter   - number of iterations (counted only when jacobian
!              reinversions are performed).
!     xi     - current mesh
!     xiold  - previous mesh
!     ipred  = 0  if relax is determined by a correction
!            = 1  if relax is determined by a prediction
!     ifreez = 0  if the jacobian is to be updated
!            = 1  if the jacobian is currently fixed (frozen)
!     iconv  = 0  if no previous convergence has been obtained
!            = 1  if convergence on a previous mesh has been obtained
!     icare  =-1  no convergence occurred (used for regular problems)
!            = 0  a regular problem
!            = 1  a sensitive problem
!            = 2  used for continuation (see description of ipar(10)
!                 in colnew).
!     rnorm  - norm of rhs (right hand side) for current iteration
!     rnold  - norm of rhs for previous iteration
!     anscl  - scaled norm of newton correction
!     anfix  - scaled norm of newton correction at next step
!     anorm  - scaled norm of a correction obtained with jacobian fixed
!     nz     - number of components of  z  (see subroutine approx)
!     ndmz   - number of components of  dmz  (see subroutine approx)
!     imesh  - a control variable for subroutines newmsh and errchk
!            = 1  the current mesh resulted from mesh selection
!                 or is the initial mesh.
!            = 2  the current mesh resulted from doubling the
!                 previous mesh
!
!**********************************************************************
!
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION XI(1), XIOLD(1), Z(1), DMZ(1), RHS(1)
      DIMENSION G(1), W(1), V(1), VALSTR(1), SLOPE(1), ACCUM(1)
      DIMENSION DELZ(1), DELDMZ(1), DQZ(1), DQDMZ(1) , FIXPNT(1)
      DIMENSION DUMMY(1), SCALE(1), DSCALE(1)
      DIMENSION INTEGS(3,1), IPVTG(1), IPVTW(1)
!
      COMMON /COLOUT/ PRECIS, IOUT, IPRINT
      COMMON /COLORD/ K, NCOMP, MSTAR, KD, MMAX, M(20)
      COMMON /COLAPR/ N, NOLD, NMAX, NZ, NDMZ
      COMMON /COLMSH/ MSHFLG, MSHNUM, MSHLMT, MSHALT
      COMMON /COLSID/ ZETA(40), ALEFT, ARIGHT, IZETA, IDUM
      COMMON /COLNLN/ NONLIN, ITER, LIMIT, ICARE, IGUESS
      COMMON /COLEST/ TOL(40), WGTMSH(40), WGTERR(40), TOLIN(40), &
                      ROOT(40), JTOL(40), LTOL(40), NTOL
!
      EXTERNAL FSUB, DFSUB, GSUB, DGSUB, GUESS
      OPTIONAL :: DFSUB, DGSUB

!---------------------------------------------------------------------------
!     Called by subroutines : Colnew
!               
!     Calls to subroutines : Newmsh (called only once in this routine) , Errchk (called only once in this routine) ,Scale, 
!                Lsyslv (call numerous times with different modes)
!
!-----------------------------------------------------------------------------------------------------------------------
!
!...  Constants for control of nonlinear iteration
!

!     relmin - minimum allowable value for relax where relax is the 
!              relaxation factor for the damped newton iteration (otherwise the
!              jacobian is considered singular).
!     rstart - initial value for relax when problem is sensitive
!     lmtfrz - maximum value for ifrz, the number of fixed jacobian iterations,  before performing a reinversion
!

      RELMIN = 1.D-3
      RSTART = 1.D-2
      LMTFRZ = 4
!
!...  Compute the maximum tolerance
!
      CHECK = 0.D0
      DO 10 I = 1, NTOL
   10   CHECK = DMAX1 ( TOLIN(I), CHECK )
!
!

!     iconv  = 0  if no previous convergence has been obtained
!     imesh  - a control variable for subroutines newmsh and errchk
!            = 1  the current mesh resulted from mesh selection
!                 or is the initial mesh.
!
      IMESH = 1
      ICONV = 0
      IF ( NONLIN .EQ. 0 ) ICONV = 1

!     What are icor=0,noconv=0 ?

      ICOR = 0
      NOCONV = 0

!     Msing = 0 implies the collocation matrix is not singular.

      MSING = 0
!
!...  The main iteration begins here .
!...  Loop 20 is executed until error tolerances are satisfied or
!...  the code fails (due to a singular matrix or storage limitations)
!     The GOTO 20 statement is at the end of this contrl routine.
!
   20      CONTINUE
!
!...       Initialization for a new mesh
!
           ! Initialize the number of iterations to 0 (counted only when Jacobian reinversions are performed)
!
           ITER = 0
           IF ( NONLIN .GT. 0 )                     GO TO 50
!
!...       The LINEAR case.
!...       Set up and solve equations
!          Call Lsyslv  with mode = 0 to set up the collocation matrices v,w,g and the right hand side and solve.
!
           CALL LSYSLV (MSING, XI, XIOLD, DUMMY, DUMMY, Z, DMZ, G, &
                W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 0, &
                FSUB, DFSUB, GSUB, DGSUB, GUESS )
!
!...       Check for a singular matrix
!          Go to line 400 and print out mesh values as we have convergence
!          and a non-singular collocation matrix.
!
           IF ( MSING .EQ. 0 )                      GO TO 400


   30      IF ( MSING .LT. 0 )                      GO TO 40

           ! A local elimination matrix is singular
           ! GOTO 460 mean pick a new mesh and try again 

           IF ( IPRINT .LT. 1 )  WRITE (IOUT,495)
           GO TO 460
   40      IF ( IPRINT .LT. 1 )  WRITE (IOUT,490)

!          The global BVP-Matrix is singular.
!          Iflag = 0 implies the collocation matrix is singular

           IFLAG = 0

!          Iflag = 1 means normal return otherwise return back to user.

           RETURN
!
!
!
!...       Iteration loop for NONLINEAR case
!...       Define the initial relaxation parameter (= relax)
!
   50      RELAX = 1.D0
!
!...       Check for previous convergence and problem sensitivity
!
!          (See icare) If a sensitive problem or no convergence occured for a regular problem reduce relax from 1 to a smaller value.

           IF ( ICARE .EQ. 1 .OR. ICARE .EQ. (-1) )  RELAX = RSTART

!          If no previous convergence has been obtained goto 160
!          where it says if:
!...       no previous convergence has been obtained. proceed
!...       with the damped newton method.
!...       evaluate rhs and find the first newton correction.

           IF ( ICONV .EQ. 0 )                      GO TO 160
!
!...       Convergence on a previous mesh has been obtained.    thus
!...       we have a very good initial approximation for the newton
!...       process. Proceed with one full newton and then iterate
!...       with a fixed jacobian.
!
!          Ifreez =0 implies jacobian is to be updated.

           IFREEZ = 0
!
!...       Evaluate right hand side and its norm  and
!...       find the first newton correction.
!          Call LSYSLV with mode = 1 implies first interation of nonlinear problems only.
!          LSYSLV solves linear systems of collocation equations.
!
           CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G, &
                W, V, RHS, DQDMZ, INTEGS, IPVTG, IPVTW, RNOLD, 1, &
                FSUB, DFSUB, GSUB, DGSUB, GUESS )
!
           IF ( IPRINT .LT. 0 )  WRITE(IOUT,530)
           IF ( IPRINT .LT. 0 )  WRITE (IOUT,510) ITER, RNOLD
           GO TO 70
!
!...       solve for the next iterate .
!...       the value of ifreez determines whether this is a full
!...       newton step (=0) or a fixed jacobian iteration (=1).
!
   60      IF ( IPRINT .LT. 0 )  WRITE (IOUT,510) ITER, RNORM
           RNOLD = RNORM


           CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G,&
                W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM,&
                3+IFREEZ, FSUB, DFSUB, GSUB, DGSUB, GUESS )
!
!...       check for a singular matrix
!
   70      IF ( MSING .NE. 0 )                      GO TO 30
           IF ( IFREEZ .EQ. 1 )                     GO TO 80
!
!...       a full newton step
!
           ITER = ITER + 1
           IFRZ = 0
   80      CONTINUE
!
!...       update   z and dmz , compute new  rhs  and its norm
!
           DO 90 I = 1, NZ
             Z(I) = Z(I) + DELZ(I)
   90      CONTINUE
           DO 100 I = 1, NDMZ
             DMZ(I) = DMZ(I) + DELDMZ(I)
  100      CONTINUE


           CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G,&
                W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 2,&
                FSUB, DFSUB, GSUB, DGSUB, GUESS )
!
!...       check monotonicity. if the norm of  rhs  gets smaller,
!...       proceed with a fixed jacobian; else proceed cautiously,
!...       as if convergence has not been obtained before (iconv=0).
!
           IF ( RNORM .LT. PRECIS )                 GO TO 390
           IF ( RNORM .GT. RNOLD )                  GO TO 130
           IF ( IFREEZ .EQ. 1 )                     GO TO 110
           IFREEZ = 1
           GO TO 60
!
!...       verify that the linear convergence with fixed jacobian
!...       is fast enough.
!
  110      IFRZ = IFRZ + 1
           IF ( IFRZ .GE. LMTFRZ )       IFREEZ = 0
           IF ( RNOLD .LT. 4.D0*RNORM )  IFREEZ = 0
!
!...       check convergence (iconv = 1).
!
           DO 120 IT = 1, NTOL
             INZ = LTOL(IT)
             DO 120 IZ = INZ, NZ, MSTAR
               IF ( DABS(DELZ(IZ)) .GT.    &
                 TOLIN(IT) * (DABS(Z(IZ)) + 1.D0))  GO TO 60
  120      CONTINUE
!
!...       convergence obtained
!

           IF ( IPRINT .LT. 1 )  WRITE (IOUT,560) ITER

           GO TO 400

!
!...      convergence of fixed jacobian iteration failed.
!
  130      IF ( IPRINT .LT. 0 )  WRITE (IOUT,510) ITER, RNORM
           IF ( IPRINT .LT. 0 )  WRITE (IOUT,540)
           ICONV = 0
           RELAX = RSTART
           DO 140 I = 1, NZ
             Z(I) = Z(I) - DELZ(I)
  140      CONTINUE
           DO 150 I = 1, NDMZ
             DMZ(I) = DMZ(I) - DELDMZ(I)
  150      CONTINUE
!
!...       update old mesh
!
           NP1 = N + 1
           DO 155 I = 1, NP1
  155        XIOLD(I) = XI(I)
           NOLD = N
!
           ITER = 0
!
!...       No previous convergence has been obtained. proceed
!...       with the damped newton method.
!...       evaluate rhs and find the first newton correction.
!
!          Print out full damped newton iteration here.

  160      IF(IPRINT .LT. 0)  WRITE (IOUT,500)


           CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G,&
                W, V, RHS, DQDMZ, INTEGS, IPVTG, IPVTW, RNOLD, 1,&
                FSUB, DFSUB, GSUB, DGSUB, GUESS )


!
!...       Check for a singular matrix. GOTO 30 could result in return to user.
!
           IF ( MSING .NE. 0 )                       GO TO 30
!
!          Matrix is not singular so do some
!...       bookkeeping for first mesh
!
           IF ( IGUESS .EQ. 1 )  IGUESS = 0
!
!...       Find initial scaling.
!          The purpose of scale is to provide a proper scaling of the state variables, used
!          to control the damping factor for a new iteration [2].
!
           CALL SKALE (N, MSTAR, KD, Z, XI, SCALE, DSCALE)
           GO TO 220
!
!...       main iteration loop
!
  170      RNOLD = RNORM
           IF ( ITER .GE. LIMIT )                   GO TO 430
!
!...       update scaling
!
           CALL SKALE (N, MSTAR, KD, Z, XI, SCALE, DSCALE)
!
!...       compute norm of newton correction with new scaling
!
           ANSCL = 0.D0
           DO 180 I = 1, NZ
             ANSCL = ANSCL + (DELZ(I) * SCALE(I))**2
  180      CONTINUE
           DO 190 I = 1, NDMZ
             ANSCL = ANSCL + (DELDMZ(I) * DSCALE(I))**2
  190      CONTINUE
           ANSCL = DSQRT(ANSCL / DFLOAT(NZ+NDMZ))
!
!...       find a newton direction


           CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G, &
                W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 3, &
                FSUB, DFSUB, GSUB, DGSUB, GUESS )
!
!...       check for a singular matrix
!
           IF ( MSING .NE. 0 )                      GO TO 30
!
!...       predict relaxation factor for newton step.
!
           ANDIF = 0.D0
           DO 200 I = 1, NZ
             ANDIF = ANDIF + ((DQZ(I) - DELZ(I)) * SCALE(I))**2
  200      CONTINUE
           DO 210 I = 1, NDMZ
             ANDIF = ANDIF + ((DQDMZ(I) - DELDMZ(I)) * DSCALE(I))**2
  210      CONTINUE
           ANDIF = DSQRT(ANDIF/DFLOAT(NZ+NDMZ) + PRECIS)
           RELAX = RELAX * ANSCL / ANDIF
           IF ( RELAX .GT. 1.D0 )  RELAX = 1.D0
  220      RLXOLD = RELAX
           IPRED = 1
           ITER = ITER + 1
!
!...       Determine a new  z and dmz  and find new  rhs  and its norm
!
           DO 230 I = 1, NZ
             Z(I) = Z(I)  +  RELAX * DELZ(I)
  230      CONTINUE
           DO 240 I = 1, NDMZ
             DMZ(I) = DMZ(I)  +  RELAX * DELDMZ(I)
  240      CONTINUE


!          Call Lsyslv with mode = 2 implies setting up the rhs only and computing its norm.


  250      CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DQZ, DQDMZ, G, &
                W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 2,&
                FSUB, DFSUB, GSUB, DGSUB, GUESS )
!
!...       compute a fixed jacobian iterate (used to control relax)
!          Call LSYSLV with mode = 4  so only forward and back ssubstition are perform.

!
           CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DQZ, DQDMZ, G, &
                W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 4,&
                FSUB, DFSUB, GSUB, DGSUB, GUESS )
!
!...       Find scaled norms of various terms used to correct Relax
!
           ANORM = 0.D0
           ANFIX = 0.D0
           DO 260 I = 1, NZ
             ANORM = ANORM  +  (DELZ(I) * SCALE(I))**2
             ANFIX = ANFIX  +  (DQZ(I) * SCALE(I))**2
  260      CONTINUE
           DO 270 I = 1, NDMZ
             ANORM = ANORM  +  (DELDMZ(I) * DSCALE(I))**2
             ANFIX = ANFIX  +  (DQDMZ(I) * DSCALE(I))**2
  270      CONTINUE
           ANORM = DSQRT(ANORM / DFLOAT(NZ+NDMZ))
           ANFIX = DSQRT(ANFIX / DFLOAT(NZ+NDMZ))
           IF ( ICOR .EQ. 1 )                         GO TO 280
           IF (IPRINT .LT. 0)  WRITE (IOUT,520) ITER, RELAX, ANORM,&
                 ANFIX, RNOLD, RNORM
           GO TO 290
  280      IF (IPRINT .LT. 0) WRITE (IOUT,550) RELAX, ANORM, ANFIX,&
                 RNOLD, RNORM
  290      ICOR = 0
!
!...       Check for monotonic decrease in  delz and deldmz.
!
           IF (ANFIX.LT.PRECIS .OR. RNORM.LT.PRECIS)  GO TO 390
           IF ( ANFIX .GT. ANORM )                    GO TO 300
!
!...       We have a decrease.
!...       If  dqz  and dqdmz  small, check for convergence
!
           IF ( ANFIX .LE. CHECK )                    GO TO 350
!
!...       Correct the predicted  relax  unless the corrected
!...       value is within 10 percent of the predicted one.
!
           IF ( IPRED .NE. 1 )                        GO TO 170
  300      IF ( ITER .GE. LIMIT )                     GO TO 430
!
!...       Correct the relaxation factor.
!
           IPRED = 0
           ARG = (ANFIX/ANORM - 1.D0) / RELAX + 1.D0
           IF ( ARG .LT. 0.D0 )                       GO TO 170
           IF (ARG .LE. .25D0*RELAX+.125D0*RELAX**2)  GO TO 310
           FACTOR = -1.D0 + DSQRT (1.D0+8.D0 * ARG)
           IF ( DABS(FACTOR-1.D0) .LT. .1D0*FACTOR )  GO TO 170
           IF ( FACTOR .LT. 0.5D0 )  FACTOR = 0.5D0
           RELAX = RELAX / FACTOR
           GO TO 320
  310      IF ( RELAX .GE. .9D0 )                     GO TO 170
           RELAX = 1.D0
  320      ICOR = 1
           IF ( RELAX .LT. RELMIN )                   GO TO 440
           FACT = RELAX - RLXOLD
           DO 330 I = 1, NZ
            Z(I) = Z(I)  +  FACT * DELZ(I)
  330      CONTINUE
           DO 340 I = 1, NDMZ
             DMZ(I) = DMZ(I)  +  FACT * DELDMZ(I)
  340      CONTINUE
           RLXOLD = RELAX
           GO TO 250
!
!...       Check convergence using tolerance values (iconv = 0).
!
  350      CONTINUE
           DO 360 IT = 1, NTOL
             INZ = LTOL(IT)
             DO 360 IZ = INZ, NZ, MSTAR
               IF ( DABS(DQZ(IZ)) .GT. &
               TOLIN(IT) * (DABS(Z(IZ)) + 1.D0) )   GO TO 170
  360      CONTINUE

!...       Convergence obtained.
!          Print out convergence after so many iterations.
!
           IF ( IPRINT .LT. 1 )  WRITE (IOUT,560) ITER
!
!...       Since convergence obtained, update  z and dmz  with term
!...       from the fixed jacobian iteration.
!
           DO 370 I = 1, NZ
             Z(I) = Z(I)  +  DQZ(I)
  370      CONTINUE
           DO 380 I = 1, NDMZ
             DMZ(I) = DMZ(I)  +  DQDMZ(I)
  380      CONTINUE


  390      IF ( (ANFIX .LT. PRECIS .OR. RNORM .LT. PRECIS) &
                 .AND. IPRINT .LT. 1 )  WRITE (IOUT,560) ITER

           ICONV = 1


           IF ( ICARE .EQ. (-1) )  ICARE = 0
!
!...       If full output has been requested, print values of the
!...       solution components   z  at the meshpoints.
!
  400      IF ( IPRINT .GE. 0 )                     GO TO 420
           DO 410 J = 1, MSTAR
             WRITE(IOUT,610) J
  410      WRITE(IOUT,620) (Z(LJ), LJ = J, NZ, MSTAR)
!
!...       Check for error tolerance satisfaction
!
  420      CONTINUE

!          CALL ERRCHK only if the current mesh resulted from doubling the old mesh.

           IF (IMESH .EQ. 2) then     
                   
                   CALL ERRCHK (XI, Z, DMZ, VALSTR, IFIN,FSUB)

           endif


           IF ( IMESH .EQ. 1 .OR. &
                IFIN .EQ. 0 .AND. ICARE .NE. 2)     GO TO 460
           IFLAG = 1
           RETURN
!
!...       Diagnostics for failure of nonlinear iteration.
!
  430      IF ( IPRINT .LT. 1 )  WRITE (IOUT,570) ITER

           GO TO 450
  440      IF( IPRINT .LT. 1 )  WRITE(IOUT,580) RELAX, RELMIN
  450      IFLAG = -2
           NOCONV = NOCONV + 1
           IF ( ICARE .EQ. 2 .AND. NOCONV .GT. 1 )  RETURN
           IF ( ICARE .EQ. 0 )  ICARE = -1
!
!...       Update old mesh
!
  460      NP1 = N + 1
           DO 470 I = 1, NP1
  470        XIOLD(I) = XI(I)
           NOLD = N
!
!...       Pick a new mesh
!...       check safeguards for mesh refinement
!
           IMESH = 1
           IF ( ICONV .EQ. 0 .OR. MSHNUM .GE. MSHLMT &
                             .OR. MSHALT .GE. MSHLMT )  IMESH = 2
           IF ( MSHALT .GE. MSHLMT .AND. &
                MSHNUM .LT. MSHLMT )  MSHALT = 1


!          Set imesh = 2 once the meshpoint error estimates satisfy the user tolerance (switch = 1) to force mesh doubling.
!          Once switch = 1, only perform mesh doubling.


           IF(K.GE.2.AND.K.LE.4.AND.SWITCH.EQ.1)THEN
              IMESH=2
           ENDIF


!           Save the values for the meshpoints, 'xiold', the solution approximations, z, and the derivative approximations, 'dmz'.
!           These will be used to compute the superconvergent solution approximations as soon as switch  = 1 in the ERRCHK routine.


            IF(K.GE.2.AND.K.LE.4)THEN


             DO I=1,NOLD+1
              MESHOLD(I)=XIOLD(I)
             ENDDO

             DO I=1,(NOLD+1)*MSTAR
              ZOLD(I)=Z(I)
             ENDDO
            
            
             DO I = 1, NDMZ
              DMZOLD(I) = DMZ(I)
             ENDDO


            ENDIF


!          Call newmsh whos purpose is to select a mesh on which a collocation solution is to be determined.

           CALL NEWMSH (IMESH, XI, XIOLD, Z, DMZ, VALSTR, &
                        SLOPE, ACCUM, NFXPNT, FIXPNT,FSUB)
!
!...       Exit if expected n is too large (but may try n=nmax once)
!
           IF ( N .LE. NMAX )                       GO TO 480
           N = N / 2
           IFLAG = -1
           IF ( ICONV .EQ. 0 .AND. IPRINT .LT. 1 )  WRITE (IOUT,590)
           IF ( ICONV .EQ. 1 .AND. IPRINT .LT. 1 )  WRITE (IOUT,600)
           RETURN
  480      IF ( ICONV .EQ. 0 )  IMESH = 1
           IF ( ICARE .EQ. 1 )  ICONV = 0

!     GOTO 20 is part of the main loop.
!...  The main iteration ends here .
!...  Loop 20 is executed until error tolerances are satisfied or
!...  the code fails (due to a singular matrix or storage limitations)
! 

      GO TO 20

!     Possibly insert conditioning constant estimation routine here.
!     This includes call to SBBLOK

!     ---------------------------------------------------------------
  490 FORMAT(//35H THE GLOBAL BVP-MATRIX IS SINGULAR )
  495 FORMAT(//40H A LOCAL ELIMINATION MATRIX IS SINGULAR )
  500 FORMAT(/30H FULL DAMPED NEWTON ITERATION,)
  510 FORMAT(13H ITERATION = , I3, 15H  NORM (RHS) = , D10.2)
  520 FORMAT(13H ITERATION = ,I3,22H  RELAXATION FACTOR = ,D10.2 &
             /33H NORM OF SCALED RHS CHANGES FROM ,D10.2,3H TO,D10.2 &
             /33H NORM   OF   RHS  CHANGES  FROM  ,D10.2,3H TO,D10.2, &
             D10.2)
  530 FORMAT(/27H FIXED JACOBIAN ITERATIONS,)
  540 FORMAT(/35H SWITCH TO DAMPED NEWTON ITERATION,)
  550 FORMAT(40H RELAXATION FACTOR CORRECTED TO RELAX = , D10.2 &
             /33H NORM OF SCALED RHS CHANGES FROM ,D10.2,3H TO,D10.2 &
             /33H NORM   OF   RHS  CHANGES  FROM  ,D10.2,3H TO,D10.2 &
             ,D10.2)
  560 FORMAT(/18H CONVERGENCE AFTER , I3,11H ITERATIONS /)
  570 FORMAT(/22H NO CONVERGENCE AFTER , I3, 11H ITERATIONS/)
  580 FORMAT(/37H NO CONVERGENCE.  RELAXATION FACTOR =,D10.3 &
             ,24H IS TOO SMALL (LESS THAN, D10.3, 1H)/)
  590 FORMAT(18H  (NO CONVERGENCE) )
  600 FORMAT(50H  (PROBABLY TOLERANCES TOO STRINGENT, OR NMAX TOO &
             ,6HSMALL) )
  610 FORMAT( 19H MESH VALUES FOR Z(, I2, 2H), )
  620 FORMAT(1H , 8D15.7)
 

     END SUBROUTINE CONTRL


      SUBROUTINE SKALE (N, MSTAR, KD, Z, XI, SCALE, DSCALE)
!
!**********************************************************************
!
!   purpose
!            provide a proper scaling of the state variables, used
!            to control the damping factor for a newton iteration [2].
!
!   variables
!
!            n      = number of mesh subintervals
!            mstar  = number of unknomns in z(u(x))
!            kd     = number of unknowns in dmz
!            z      = the global unknown vector
!            xi     = the current mesh
!            scale  = scaling vector for z
!            dscale = scaling vector for dmz
!
!**********************************************************************
!
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION Z(MSTAR,1), SCALE(MSTAR,1), DSCALE(KD,1)
      DIMENSION XI(1), BASM(5)
!
      COMMON /COLORD/ K, NCOMP, ID1, ID2, MMAX, M(20)

!-----------------------------------------------------------------------------------------------------------------------
!     Called by subroutines : contrl
!               
!     Calls to : No subroutine calls 
!
!-----------------------------------------------------------------------------------------------------------------------
!   
!
      BASM(1) = 1.D0
      DO 50 J=1,N
        IZ = 1
        H = XI(J+1) - XI(J)
        DO 10 L = 1, MMAX
          BASM(L+1) = BASM(L) * H / DFLOAT(L)
  10    CONTINUE
        DO 40 ICOMP = 1, NCOMP
          SCAL = (DABS(Z(IZ,J)) + DABS(Z(IZ,J+1))) * .5D0 + 1.D0
          MJ = M(ICOMP)
          DO 20 L = 1, MJ
            SCALE(IZ,J) = BASM(L) / SCAL
            IZ = IZ + 1
  20      CONTINUE
          SCAL = BASM(MJ+1) / SCAL
          DO 30 IDMZ = ICOMP, KD, NCOMP
            DSCALE(IDMZ,J) = SCAL
  30      CONTINUE
  40    CONTINUE
  50  CONTINUE
      NP1 = N + 1
      DO 60 IZ = 1, MSTAR
        SCALE(IZ,NP1) = SCALE(IZ,N)
  60  CONTINUE
      RETURN
 

     END SUBROUTINE SKALE




!----------------------------------------------------------------------
!                            p a r t  2
!          mesh selection, error estimation, (and related
!          constant assignment) routines -- see [3], [4]
!----------------------------------------------------------------------
!
      SUBROUTINE NEWMSH (MODE, XI, XIOLD, Z, DMZ, VALSTR, &
                         SLOPE, ACCUM, NFXPNT, FIXPNT,FSUB)
!
!**********************************************************************
!
!   purpose
!            select a mesh on which a collocation solution is to be
!            determined
!
!  
!             Nfixpnt-often 0 extra points given inside the boundary conditions
!             Fixpnt - Often none
!             Fsub - Used supplied function routine
!
!
!                           there are 5 possible modes of action:
!            mode = 5,4,3 - deal mainly with definition of an initial
!                           mesh for the current boundary value problem
!                 = 2,1   - deal with definition of a new mesh, either
!                           by simple mesh halving or by mesh selection
!                           (imesh=imode = either 1,2 when newmsh is being called by contrl.
!            more specifically, for
!            mode = 5  an initial (generally nonuniform) mesh is
!                      defined by the user and no mesh selection is to
!                      be performed
!                 = 4  an initial (generally nonuniform) mesh is
!                      defined by the user
!                 = 3  a simple uniform mesh (except possibly for some
!                      fixed points) is defined; n= no. of subintervals
!                 = 1  the automatic mesh selection procedure is used
!                      (see [3] for details)
!                 = 2  a simple mesh halving is performed
!
!    
!**********************************************************************
!
!   variables
!
!            n      = number of mesh subintervals
!            nold   = number of subintervals for former mesh
!            xi     - mesh point array
!            xiold  - former mesh point array
!            mshlmt - maximum no. of mesh selections which are permitted
!                     for a given n before mesh halving
!            mshnum - no. of mesh selections which have actually been
!                     performed for the given n
!            mshalt - no. of consecutive times ( plus 1 ) the mesh
!                     selection has alternately halved and doubled n.
!                     if mshalt .ge. mshlmt then  contrl  requires
!                     that the current mesh be halved.
!            mshflg = 1  the mesh is a halving of its former mesh
!                       (so an error estimate has been calculated)
!                   = 0  otherwise
!            iguess - ipar(9) in subroutine colnew.  it is used
!                     here only for mode=5 and 4, where
!                   = 2 the subroutine sets xi=xiold.  this is
!                       used e.g. if continuation is being per-
!                       formed, and a mesh for the old differen-
!                       tial equation is being used
!                   = 3 same as for =2, except xi uses every other
!                       point of xiold (so mesh xiold is mesh xi
!                       halved)
!                   = 4 xi has been defined by the user, and an old
!                       mesh xiold is also available
!                       otherwise, xi has been defined by the user
!                       and we set xiold=xi in this subroutine
!            slope  - an approximate quantity to be equidistributed for
!                     mesh selection (see [3]), viz,
!                             .                        (k+mj)
!                     slope(i)=     max   (weight(l) *u      (xi(i)))
!                               1.le.l.le.ntol         j
!
!                     where j=jtol(l)
!            slphmx - maximum of slope(i)*(xiold(i+1)-xiold(i)) for
!                     i = 1 ,..., nold.
!            accum  - accum(i) is the integral of  slope  from  aleft
!                     to  xiold(i).
!            valstr - is assigned values needed in  errchk  for the
!                     error estimate.
!**********************************************************************
!
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION D1(40), D2(40), SLOPE(1), ACCUM(*), VALSTR(1)
      DIMENSION XI(1), XIOLD(*), Z(1), DMZ(1), FIXPNT(1), DUMMY(1)
!
!     'y_super' - yi values used for the SCSA
!
!     'sk_discrete_super' - discrete stages for the SCSA
!
!     'sk_interp_super' - extra interpolant stages for SCSA
!

      DOUBLE PRECISION Y_SUPER(160*(3000+1))
      DOUBLE PRECISION SK_DISCRETE_SUPER(3000*10*160)
      DOUBLE PRECISION SK_INTERP_SUPER(3000*10*160)

!     's_super' - the number of discrete stages in the SCSA
!
!     's_star_super' - the number of extra interpolant stages in the SCSA
!

      INTEGER S_SUPER
      INTEGER S_STAR_SUPER

      COMMON /STAGES/ Y_SUPER, SK_DISCRETE_SUPER, SK_INTERP_SUPER

      COMMON /COLOUT/ PRECIS, IOUT, IPRINT
      COMMON /COLORD/ K, NCOMP, MSTAR, KD, MMAX, M(20)
      COMMON /COLAPR/ N, NOLD, NMAX, NZ, NDMZ
      COMMON /COLMSH/ MSHFLG, MSHNUM, MSHLMT, MSHALT
      COMMON /COLNLN/ NONLIN, ITER, LIMIT, ICARE, IGUESS
      COMMON /COLSID/  ZETA(40), ALEFT, ARIGHT, IZETA, IDUM
      COMMON /COLBAS/ B(28), ACOL(28,7), ASAVE(28,7)
      COMMON /COLEST/ TOL(40), WGTMSH(40), WGTERR(40), TOLIN(40), &
                      ROOT(40), JTOL(40), LTOL(40), NTOL

      EXTERNAL FSUB

!----------------------------------------------------------------------------------------------------------------------
!
!     Called by subroutines: colnew, contrl
!               
!     Calls to subroutines : horder, approx
!
!-----------------------------------------------------------------------------------------------------------------------
!
!
!     The correct values for 's_super' and 's_star_super' for the SCSA depend on the value of 'k', the number of collocation points 

      IF(K.EQ.2)THEN
        S_SUPER=4
        S_STAR_SUPER=4
      ENDIF
      
      IF(K.EQ.3)THEN
        S_SUPER=5
        S_STAR_SUPER=6
      ENDIF

      IF(K.EQ.4)THEN
        S_SUPER=6
        S_STAR_SUPER=9
      ENDIF

!
!
      NFXP1 = NFXPNT +1

!     For first time called by newmsh mode is likely 3 for generating a uniform mesh.
!     So GOTO line 50.





      GO TO (180, 100, 50, 20, 10), MODE
!
!...  mode=5   set mshlmt=1 so that no mesh selection is performed
!
   10 MSHLMT = 1
!
!...  mode=4   the user-specified initial mesh is already in place.
!
   20 IF ( IGUESS .LT. 2 )                          GO TO 40
!
!...  iguess=2, 3 or 4.
!
      NOLDP1 = NOLD + 1
      IF (IPRINT .LT. 1)  WRITE(IOUT,360) NOLD, (XIOLD(I), I=1,NOLDP1)
      IF ( IGUESS .NE. 3 )                          GO TO 40
!
!...  if iread ( ipar(8) ) .ge. 1 and iguess ( ipar(9) ) .eq. 3
!...  then the first mesh is every second point of the
!...  mesh in  xiold .
!
      N = NOLD /2
      I = 0
      DO 30 J = 1, NOLD, 2
           I = I + 1
   30 XI(I) = XIOLD(J)
   40 CONTINUE
      NP1 = N + 1
      XI(1) = ALEFT
      XI(NP1) = ARIGHT
      GO TO 320
!
!...  Mode=3   generate a (piecewise) uniform mesh. if there are
!...  fixed points then ensure that the n being used is large enough.
!     Loop from 50 to 90 then go to end of newmsh using GOTO 320.
!
   50 IF ( N .LT. NFXP1 )  N = NFXP1
      NP1 = N + 1
      XI(1) = ALEFT
      ILEFT = 1
      XLEFT = ALEFT
!
!...  Loop over the subregions between fixed points.
!
      DO 90 J = 1, NFXP1
           XRIGHT = ARIGHT
           IRIGHT = NP1
           IF ( J .EQ. NFXP1 )                      GO TO 60
           XRIGHT = FIXPNT(J)
!
!...       Determine where the j-th fixed point should fall in the
!...       new mesh - this is xi(iright) and the (j-1)st fixed
!...       point is in xi(ileft)
!
           NMIN = (XRIGHT-ALEFT) / (ARIGHT-ALEFT) * DFLOAT(N) + 1.5D0
           IF (NMIN .GT. N-NFXPNT+J)  NMIN = N - NFXPNT + J
           IRIGHT = MAX0 (ILEFT+1, NMIN)
   60      XI(IRIGHT) = XRIGHT
!
!...       Generate equally spaced points between the j-1st and the
!...       j-th fixed points.
!
           NREGN = IRIGHT - ILEFT - 1
           IF ( NREGN .EQ. 0 )                      GO TO 80
           DX = (XRIGHT - XLEFT) / DFLOAT(NREGN+1)
           DO 70 I = 1, NREGN
   70      XI(ILEFT+I) = XLEFT  +  DFLOAT(I) * DX
   80      ILEFT = IRIGHT
           XLEFT = XRIGHT
   90 CONTINUE
      GO TO 320
!
!...  Mode=2  halve the current mesh (i.e. double its size)
!
  100 N2 = 2 * N
!
!...  Check that n does not exceed storage limitations
!

      IF ( N2 .LE. NMAX )                           GO TO 120
!
!...  If possible, try with n=nmax. redistribute first.
!
      IF ( MODE .EQ. 2 )                            GO TO 110
      N = NMAX / 2
      GO TO 220
  110 IF ( IPRINT .LT. 1 )  WRITE(IOUT,370)
      N = N2
      RETURN

!...  Calculate the old approximate solution values at
!...  points to be used in  errchk  for error estimates.
!...  if  mshflg  = 1 an error estimate was obtained for
!...  for the old approximation so half the needed values
!...  will already be in  valstr .
!

  120    CONTINUE
  
!
      IF ( MSHFLG .EQ. 0)                          GO TO 140
!
!...  Save in  valstr  the values of the old solution
!...  at the relative positions 1/6 and 5/6 in each subinterval.
!
!     The 'evaluate' subroutine is used to compute the SCSAs at 1/3 and 2/3 across the subinterval 
!     for the primary error estimates in the ERRCHK routine.
!
!     This is necessary for 'k' = 2,3 or 4 during the last part of the computation ('switch' = 1 ) 
!
      IF(K.GE.2.AND.K.LE.4.AND.SWITCH.EQ.1)THEN
      

        KSTORE = 1

        DO 130 I = 1, NOLD
         HD6 = (XIOLD(I+1) - XIOLD(I)) / 6.D0
         X = XIOLD(I) + HD6

         CALL EVALUATE(MSTAR,NOLD,XIOLD,Y_SUPER,X,S_SUPER,S_STAR_SUPER,&
                       VALSTR(KSTORE),SK_DISCRETE_SUPER,SK_INTERP_SUPER)

         X = X + 4.D0 * HD6
         KSTORE = KSTORE  +  3 * MSTAR 

       
         CALL EVALUATE(MSTAR,NOLD,XIOLD,Y_SUPER,X,S_SUPER,S_STAR_SUPER,&
                       VALSTR(KSTORE),SK_DISCRETE_SUPER,SK_INTERP_SUPER)

         KSTORE = KSTORE  +  MSTAR
  130   CONTINUE
        GO TO 160
!
!...  Save in  valstr  the values of the old solution
!...  at the relative positions 1/6, 2/6, 4/6 and 5/6 in
!...  each subinterval.

!     The 'approx' subroutine is used to compute the primary meshpoint error estimates in the ERRCHK routine.
!
!     This is necessary for 'k' = 2,3 or 4 during the first part of the computation ('switch' = 0 ) 
!
      ELSE IF (K.GE.2.AND.K.LE.4.AND.SWITCH.EQ.0) THEN

        YKSTORE = 1
        
        DO 131 I = 1, NOLD

        X= XIOLD(I)
      
        CALL APPROX(I,X,Y_MESHOLD(YKSTORE),ASAVE(1,5),DUMMY,XIOLD,&
            NOLD,Z,DMZ,K,NCOMP,MMAX,M,MSTAR,4,DUMMY,0)       

        YKSTORE = YKSTORE  +  MSTAR

        IF(I.EQ.NOLD)THEN

            X=XIOLD(I+1)

            CALL APPROX(I,X,Y_MESHOLD(YKSTORE),ASAVE(1,7),DUMMY,XIOLD,&
            NOLD,Z,DMZ,K,NCOMP,MMAX,M,MSTAR,4,DUMMY,0)       

        ENDIF


  131   CONTINUE


        GO TO 160
!
!...  Save in  valstr  the values of the old solution
!...  at the relative positions 1/6, 2/6, 4/6 and 5/6 in
!...  each subinterval.
!
!     For 'k' = 1,5,6 and 7, use the original colnew code
!
      ELSE ! K=1,5,6 OR 7


          KSTORE = 1

          DO 132 I = 1, NOLD
           HD6 = (XIOLD(I+1) - XIOLD(I)) / 6.D0
           X = XIOLD(I) + HD6
           CALL APPROX (I, X, VALSTR(KSTORE), ASAVE(1,1), DUMMY, XIOLD, &
              NOLD, Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 4, DUMMY, 0)
           X = X + 4.D0 * HD6
           KSTORE = KSTORE  +  3 * MSTAR
           CALL APPROX (I, X, VALSTR(KSTORE), ASAVE(1,4), DUMMY, XIOLD, &
              NOLD, Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 4, DUMMY, 0)
           KSTORE = KSTORE  +  MSTAR
  132     CONTINUE

          GO TO 160    

      ENDIF 
!
  140 CONTINUE


      IF(K.GE.2.AND.K.LE.4.AND.SWITCH.EQ.1)THEN


        CALL SETUP_PHASE1(NOLD,XIOLD,ASAVE(1,3),Z,DMZ,K,NCOMP,MMAX, &
              M,MSTAR,Y_SUPER,SK_DISCRETE_SUPER,SK_INTERP_SUPER,FSUB)

        KSTORE=1

        DO 150 I = 1, N
         X = XI(I)
         HD6 = (XI(I+1) - XI(I)) / 6.D0

         DO 150 J = 1, 4
           X = X + HD6
           IF ( J.EQ.3 )  X = X + HD6 

           CALL EVALUATE(MSTAR,NOLD,XIOLD,Y_SUPER,X,S_SUPER,S_STAR_SUPER,&
                VALSTR(KSTORE),SK_DISCRETE_SUPER,SK_INTERP_SUPER)

           KSTORE = KSTORE  +  MSTAR
  150   CONTINUE


       ELSE IF (K.GE.2.AND.K.LE.4.AND.SWITCH.EQ.0) THEN
 
         YKSTORE=1

       
         DO 151 I = 1, N

           X = XI(I)

           CALL APPROX(I,X,Y_MESHOLD(YKSTORE),ASAVE(1,5),DUMMY,XIOLD,&
            NOLD,Z,DMZ,K,NCOMP,MMAX,M,MSTAR,4,DUMMY,0)       
 
            YKSTORE = YKSTORE  +  MSTAR

            IF(I.EQ.N)THEN

            X=XI(I+1)

            CALL APPROX(I,X,Y_MESHOLD(YKSTORE),ASAVE(1,7),DUMMY,XIOLD,&
            NOLD,Z,DMZ,K,NCOMP,MMAX,M,MSTAR,4,DUMMY,0)       


            ENDIF


  151    CONTINUE




       ELSE !  'K' = 1,5,6,7

           KSTORE = 1

           DO 152 I = 1, N

            X = XI(I)
            HD6 = (XI(I+1) - XI(I)) / 6.D0

            DO 152 J = 1, 4

             X = X + HD6
             IF ( J.EQ.3 )  X = X + HD6

             CALL APPROX (I, X, VALSTR(KSTORE), ASAVE(1,J), DUMMY, XIOLD, &
                           NOLD, Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 4, DUMMY, 0)
             KSTORE = KSTORE  +  MSTAR
  152     CONTINUE      




      ENDIF


  160 MSHFLG = 0
      MSHNUM = 1
      MODE = 2
!
!...  Generate the halved mesh.
!
      J = 2
      DO 170 I = 1, N
           XI(J) = (XIOLD(I) + XIOLD(I+1)) / 2.D0
           XI(J+1) = XIOLD(I+1)
  170 J = J + 2
      N = N2
      GO TO 320
!
!...  Mode=1  we do mesh selection if it is deemed worthwhile
!
  180 IF ( NOLD .EQ. 1 )                            GO TO 100
      IF ( NOLD .LE. 2*NFXPNT )                     GO TO 100
!
!...  The first interval has to be treated separately from the
!...  other intervals (generally the solution on the (i-1)st and ith
!...  intervals will be used to approximate the needed derivative, but
!...  here the 1st and second intervals are used.)
!
      I = 1
      HIOLD = XIOLD(2) - XIOLD(1)

!     Call horder to determine the highest order (piecewise constant) derivatives of the
!     current solution on the first subinterval.

      CALL HORDER (1, D1, HIOLD, DMZ, NCOMP, K)

      HIOLD = XIOLD(3) - XIOLD(2)

!     Call horder to determine the highest order (piecewise constant) derivatives of the
!     current solution on the second subinterval.

      CALL HORDER (2, D2, HIOLD, DMZ, NCOMP, K)
      ACCUM(1) = 0.D0
      SLOPE(1) = 0.D0
      ONEOVH = 2.D0 / ( XIOLD(3) - XIOLD(1) )
      DO 190 J = 1, NTOL
           JJ = JTOL(J)
           JZ = LTOL(J)
  190 SLOPE(1) = DMAX1(SLOPE(1),(DABS(D2(JJ)-D1(JJ))*WGTMSH(J)* &
                 ONEOVH / (1.D0 + DABS(Z(JZ)))) **ROOT(J))
      SLPHMX = SLOPE(1) * (XIOLD(2) - XIOLD(1))
      ACCUM(2) = SLPHMX
      IFLIP = 1
!
!...  Go through the remaining intervals generating  slope
!...  and  accum .
!
      DO 210 I = 2, NOLD
           HIOLD = XIOLD(I+1) - XIOLD(I)
           IF ( IFLIP .EQ. -1 ) &
                         CALL HORDER ( I, D1, HIOLD, DMZ, NCOMP, K)
           IF ( IFLIP .EQ. 1 ) &
                         CALL HORDER ( I, D2, HIOLD, DMZ, NCOMP, K)
           ONEOVH = 2.D0 / ( XIOLD(I+1) - XIOLD(I-1) )
           SLOPE(I) = 0.D0
!
!...       Evaluate function to be equidistributed
!
           DO 200 J = 1, NTOL
             JJ = JTOL(J)
             JZ = LTOL(J)  +  (I-1) * MSTAR
  200      SLOPE(I) = DMAX1(SLOPE(I),(DABS(D2(JJ)-D1(JJ))*WGTMSH(J)* &
                      ONEOVH / (1.D0 + DABS(Z(JZ)))) **ROOT(J))
!
!...       Accumulate approximate integral of function to be
!...       equidistributed
!
           TEMP = SLOPE(I) * (XIOLD(I+1)-XIOLD(I))
           SLPHMX = DMAX1(SLPHMX,TEMP)
           ACCUM(I+1) = ACCUM(I) + TEMP
  210 IFLIP = - IFLIP
      AVRG = ACCUM(NOLD+1) / DFLOAT(NOLD)
      DEGEQU = AVRG / DMAX1(SLPHMX,PRECIS)
!
!...  Naccum=expected n to achieve .1x user requested tolerances
!
      NACCUM = ACCUM(NOLD+1) + 1.D0

      IF ( IPRINT .LT. 0 )  WRITE(IOUT,350) DEGEQU, NACCUM

!...  Decide if mesh selection is worthwhile (otherwise, halve)
!
!     Go to 100 means halve the current mesh.


      IF ( AVRG .LT. PRECIS ) GO TO 100
      IF ( DEGEQU .GE. .5D0 ) GO TO 100

!
!...  Nmx assures mesh has at least half as many subintervals as the
!...  previous mesh
!
      NMX = MAX0 ( NOLD+1, NACCUM ) / 2
!
!...  This assures that halving will be possible later (for error est)
!
      NMAX2 = NMAX / 2
!
!...  The mesh is at most halved
!
      N = MIN0 ( NMAX2, NOLD, NMX )
  220 NOLDP1 = NOLD + 1
      IF ( N .LT. NFXP1 )  N = NFXP1
      MSHNUM = MSHNUM + 1
!
!...  If the new mesh is smaller than the old mesh set mshnum
!...  so that the next call to  newmsh  will produce a halved
!...  mesh. if n .eq. nold / 2 increment mshalt so there can not
!...  be an infinite loop alternating between n and n/2 points.
!
      IF ( N .LT. NOLD )  MSHNUM = MSHLMT
      IF ( N .GT. NOLD/2 )  MSHALT = 1
      IF ( N .EQ. NOLD/2 )  MSHALT = MSHALT + 1
      MSHFLG = 0
!
!...  Having decided to generate a new mesh with n subintervals we now
!...  do so, taking into account that the nfxpnt points in the array
!...  fixpnt must be included in the new mesh.
!
      IN = 1
      ACCL = 0.D0
      LOLD = 2
      XI(1) = ALEFT
      XI(N+1) = ARIGHT
      DO 310 I = 1, NFXP1
           IF ( I .EQ. NFXP1 )                      GO TO 250
           DO 230 J = LOLD, NOLDP1
             LNEW = J
             IF ( FIXPNT(I) .LE. XIOLD(J) )         GO TO 240
  230      CONTINUE
  240      CONTINUE
           ACCR = ACCUM(LNEW) + (FIXPNT(I)-XIOLD(LNEW))*SLOPE(LNEW-1)
           NREGN = (ACCR-ACCL) / ACCUM(NOLDP1) * DFLOAT(N) - .5D0
           NREGN = MIN0(NREGN, N - IN - NFXP1 + I)
           XI(IN + NREGN + 1) = FIXPNT(I)
           GO TO 260
  250      ACCR = ACCUM(NOLDP1)
           LNEW = NOLDP1
           NREGN = N - IN
  260      IF ( NREGN .EQ. 0 )                      GO TO 300
           TEMP = ACCL
           TSUM = (ACCR - ACCL) / DFLOAT(NREGN+1)
           DO 290 J = 1, NREGN
             IN = IN + 1
             TEMP = TEMP + TSUM
             DO 270 L = LOLD, LNEW
               LCARRY = L
               IF ( TEMP .LE. ACCUM(L) )            GO TO 280
  270        CONTINUE
  280        CONTINUE
             LOLD = LCARRY
  290      XI(IN) = XIOLD(LOLD-1) + (TEMP - ACCUM(LOLD-1)) / &
           SLOPE(LOLD-1)
  300      IN = IN + 1
           ACCL = ACCR
           LOLD = LNEW
  310 CONTINUE
      MODE = 1
  320 CONTINUE
      NP1 = N + 1

      IF ( IPRINT .LT. 1 )  WRITE(IOUT,340) N, (XI(I),I=1,NP1)

      NZ   = MSTAR * (N + 1)
      NDMZ = KD * N
      RETURN
!----------------------------------------------------------------
  340 FORMAT(/17H THE NEW MESH (OF,I5,16H SUBINTERVALS), ,100(/8F12.6))
  350 FORMAT(/21H MESH SELECTION INFO,/30H DEGREE OF EQUIDISTRIBUTION = &
             , F8.5, 28H PREDICTION FOR REQUIRED N = , I8)
  360 FORMAT(/20H THE FORMER MESH (OF,I5,15H SUBINTERVALS),, &
             100(/8F12.6))
  370 FORMAT (/23H  EXPECTED N TOO LARGE  )
 

     END SUBROUTINE NEWMSH


      SUBROUTINE CONSTS (K, RHO, COEF)
!
!**********************************************************************
!
!   purpose
!            assign (once) values to various array constants.
!
!   arrays assigned during compilation:
!     cnsts1 - weights for extrapolation error estimate
!     cnsts2 - weights for mesh selection
!              (the above weights come from the theoretical form for
!              the collocation error -- see [3])
!
!   arrays assigned during execution:
!     wgterr - the particular values of cnsts1 used for current run
!              (depending on k, m)
!     wgtmsh - gotten from the values of cnsts2 which in turn are
!              the constants in the theoretical expression for the
!              errors. the quantities in wgtmsh are 10x the values
!              in cnsts2 so that the mesh selection algorithm
!              is aiming for errors .1x as large as the user
!              requested tolerances.
!     jtol   - components of differential system to which tolerances
!              refer (viz, if ltol(i) refers to a derivative of u(j),
!              then jtol(i)=j)
!     root   - reciprocals of expected rates of convergence of compo-
!              nents of z(j) for which tolerances are specified
!     rho    - the k collocation points on (0,1)
!     coef   - are the precomputed derivatives of the basis
!     acol  -  the runge-kutta coefficients values at collocation
!              points
!
!**********************************************************************
!
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION RHO(7), COEF(K,1), CNSTS1(28), CNSTS2(28), DUMMY(1)
!
      COMMON /COLORD/ KDUM, NCOMP, MSTAR, KD, MMAX, M(20)
      COMMON /COLBAS/ B(28), ACOL(28,7), ASAVE(28,7)
      COMMON /COLEST/ TOL(40), WGTMSH(40), WGTERR(40), TOLIN(40), &
                      ROOT(40), JTOL(40), LTOL(40), NTOL
!----------------------------------------------------------------------------------------------------------------------
!
!     Called by subroutines: colnew
!               
!     Calls to subroutines: vmonde, rkbas
!
!-----------------------------------------------------------------------------------------------------------------------

!    cnsts1 - weights for extrapolation error estimate
!    cnsts2 - weights for mesh selection
!             (the above weights come from the theoretical form for
!              the collocation error -- see [3])
!
!   Both CNSTS1 and CNSTS2 have 28 values


      DATA CNSTS1 /    .25D0,     .625D-1,  7.2169D-2, 1.8342D-2,&
           1.9065D-2, 5.8190D-2, 5.4658D-3, 5.3370D-3, 1.8890D-2,&
           2.7792D-2, 1.6095D-3, 1.4964D-3, 7.5938D-3, 5.7573D-3,&
           1.8342D-2, 4.673D-3,  4.150D-4,  1.919D-3,  1.468D-3,&
           6.371D-3,  4.610D-3,  1.342D-4,  1.138D-4,  4.889D-4,&
           4.177D-4,  1.374D-3,  1.654D-3,  2.863D-3  /
      DATA CNSTS2 /   1.25D-1,   2.604D-3,  8.019D-3,  2.170D-5,&
           7.453D-5,  5.208D-4,  9.689D-8,  3.689D-7,  3.100D-6,&
           2.451D-5,  2.691D-10, 1.120D-9,  1.076D-8,  9.405D-8,&
           1.033D-6,  5.097D-13, 2.290D-12, 2.446D-11, 2.331D-10,&
           2.936D-9,  3.593D-8,  7.001D-16, 3.363D-15, 3.921D-14,&
           4.028D-13, 5.646D-12, 7.531D-11, 1.129D-9  /
!
!...  Assign weights for error estimate
!
!     
!     Wgterr is assigned- the particular values of cnsts1 used for current run
!     (depending on k, m)
!






      IF(K.EQ.2)THEN
         CNSTS1(3)=CNSTS1(2)
      ENDIF
      IF(K.EQ.3)THEN
         CNSTS1(6)=CNSTS1(4)
      ENDIF
      IF(K.EQ.4)THEN
         CNSTS1(10)=CNSTS1(7)
      ENDIF







      KOFF = K * ( K + 1 ) / 2
      IZ = 1
      DO 10 J = 1, NCOMP
           MJ = M(J)
           DO 10 L = 1, MJ
             WGTERR(IZ) = CNSTS1(KOFF - MJ + L)
             IZ = IZ + 1
   10 CONTINUE
!
!...  Assign array values for mesh selection: wgtmsh, jtol, and root
!
!
!     Wgtmsh is gotten from the values of cnsts2 which in turn are
!              the constants in the theoretical expression for the
!              errors. the quantities in wgtmsh are 10x the values
!              in cnsts2 so that the mesh selection algorithm
!              is aiming for errors .1x as large as the user
!              requested tolerances.

!     Root   - reciprocals of expected rates of convergence of compo-
!                         nents of z(j) for which tolerances are specified
!
!     Jtol   - components of differential system to which tolerances
!
      JCOMP = 1
      MTOT = M(1)
      DO 40 I = 1, NTOL
           LTOLI = LTOL(I)
   20      CONTINUE
           IF ( LTOLI .LE. MTOT )                   GO TO 30
           JCOMP = JCOMP + 1
           MTOT = MTOT + M(JCOMP)
           GO TO 20
   30      CONTINUE
           JTOL(I) = JCOMP
           WGTMSH(I) = 1.D1 * CNSTS2(KOFF+LTOLI-MTOT) / TOLIN(I)




         

          IF(K.GE.2.AND.K.LE.4)THEN

              ROOT(I)=1.0D0/( DFLOAT(2*K) )

          ELSE

             ROOT(I) = 1.D0 / DFLOAT(K+MTOT-LTOLI+1)
        
          ENDIF 





   40 CONTINUE
!
!...  Specify collocation points dependent on value of K in the following GOTO statement.
!
      GO TO (50,60,70,80,90,100,110), K
   50 RHO(1) = 0.D0
      GO TO 120
   60 RHO(2) = .57735026918962576451D0
      RHO(1) = - RHO(2)
      GO TO 120
   70 RHO(3) = .77459666924148337704D0
      RHO(2) = .0D0
      RHO(1) = - RHO(3)
      GO TO 120
   80 RHO(4) = .86113631159405257523D0
      RHO(3) = .33998104358485626480D0
      RHO(2) = - RHO(3)
      RHO(1) = - RHO(4)
      GO TO 120
   90 RHO(5) = .90617984593866399280D0
      RHO(4) = .53846931010568309104D0
      RHO(3) = .0D0
      RHO(2) = - RHO(4)
      RHO(1) = - RHO(5)
      GO TO 120
  100 RHO(6) = .93246951420315202781D0
      RHO(5) = .66120938646626451366D0
      RHO(4) = .23861918608319690863D0
      RHO(3) = -RHO(4)
      RHO(2) = -RHO(5)
      RHO(1) = -RHO(6)
      GO TO 120
  110 RHO(7) = .949107991234275852452D0
      RHO(6) = .74153118559939443986D0
      RHO(5) = .40584515137739716690D0
      RHO(4) = 0.D0
      RHO(3) = -RHO(5)
      RHO(2) = -RHO(6)
      RHO(1) = -RHO(7)
  120 CONTINUE
!
!     Some of the values for RHO are from (-1,0)
!...  map (-1,1) to (0,1) by  t = .5 * (1. + x)
!
      DO 130 J = 1, K
         RHO(J) = .5D0 * (1.D0 + RHO(J))
  130 CONTINUE
!
!...  Now find runge-kutta coeffitients b, acol and asave
!...  the values of asave are to be used in  newmsh  and errchk .
!
      DO 140 J = 1, K
         DO 135 I = 1, K
  135      COEF(I,J) = 0.D0
         COEF(J,J) = 1.D0

!     Vmonde computes the coef value given Rho and K as input.

         CALL VMONDE (RHO, COEF(1,J), K)
  140 CONTINUE

!
!   RKBAS will evaluate mesh independent runge-kutta basis for given s
!
!     s      - is the first argument, i.e. the relative position for which
!              the basis is to be evaluated ( 0. .le. s .le. 1. ).
!     coef   is inputted as  precomputed derivatives of the basis
!     k      is inputted as number of collocatin points per subinterval
!     m      is inputted as maximal order of the differential equation
!     rkb    is outputted as the the runge-kutta basis (0-th to (m-1)-th derivatives )
!     dm     is not computed here (dummy) as basis elements for m-th derivative because mode =0.

      CALL RKBAS ( 1.D0, COEF, K, MMAX, B, DUMMY, 0)
      DO 150 I = 1, K
         CALL RKBAS ( RHO(I), COEF, K, MMAX, ACOL(1,I), DUMMY, 0)
  150 CONTINUE

      CALL RKBAS ( 1.D0/6.D0, COEF, K, MMAX, ASAVE(1,1), DUMMY, 0)
      CALL RKBAS ( 1.D0/3.D0, COEF, K, MMAX, ASAVE(1,2), DUMMY, 0)
      CALL RKBAS ( 2.D0/3.D0, COEF, K, MMAX, ASAVE(1,3), DUMMY, 0)
      CALL RKBAS ( 5.D0/6.D0, COEF, K, MMAX, ASAVE(1,4), DUMMY, 0)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      CALL RKBAS (0.0d0, COEF, K, MMAX, ASAVE(1,5), DUMMY, 0)
      CALL RKBAS (1.0d0/2.0d0, COEF, K, MMAX, ASAVE(1,6), DUMMY, 0)
      CALL RKBAS (1.0d0, COEF, K, MMAX, ASAVE(1,7), DUMMY, 0)



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


      RETURN


      END SUBROUTINE CONSTS



      SUBROUTINE ERRCHK (XI, Z, DMZ, VALSTR, IFIN,FSUB)
!
!**********************************************************************
!
!      purpose
!               determine the error estimates and test to see if the
!               error tolerances are satisfied.
!
!      variables
!        xi     - current mesh points
!        valstr - values of the previous solution which are needed
!                 for the extrapolation- like error estimate.
!        wgterr - weights used in the extrapolation-like error
!                 estimate. the array values are assigned in
!                 subroutine  consts.
!        errest - storage for error estimates
!        err    - temporary storage used for error estimates
!        z      - approximate solution on mesh xi
!        ifin   - a 0-1 variable. on return it indicates whether
!                 the error tolerances were satisfied
!        mshflg - is set by errchk to indicate to newmsh whether
!                 any values of the current solution are stored in
!                 the array valstr. (0 for no, 1 for yes)
!
!**********************************************************************
!
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION ERR(40), ERREST(40), DUMMY(1)
      DIMENSION ERR1(40),ERR3(40)
      DIMENSION XI(1), Z(1), DMZ(1), VALSTR(1)

      DOUBLE PRECISION Y_SUPER(160*(3000+1))
      DOUBLE PRECISION SK_DISCRETE_SUPER(3000*10*160)
      DOUBLE PRECISION SK_INTERP_SUPER(3000*10*160)

      INTEGER S_SUPER
      INTEGER S_STAR_SUPER


      COMMON /STAGES/ Y_SUPER, SK_DISCRETE_SUPER, SK_INTERP_SUPER


      COMMON /COLOUT/ PRECIS, IOUT, IPRINT
      COMMON /COLORD/ K, NCOMP, MSTAR, KD, MMAX, M(20)
      COMMON /COLAPR/ N, NOLD, NMAX, NZ, NDMZ
      COMMON /COLMSH/ MSHFLG, MSHNUM, MSHLMT, MSHALT
      COMMON /COLBAS/ B(28), ACOL(28,7), ASAVE(28,7)
      COMMON /COLEST/ TOL(40), WGTMSH(40), WGTERR(40), TOLIN(40), &
                      ROOT(40), JTOL(40), LTOL(40), NTOL


      EXTERNAL FSUB

!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Contrl
!               
!     Calls to : Approx
!
!-----------------------------------------------------------------------------------------------------------------------



!     
!    Using meshpoint error control for  2 <= k <= 4 for first part of computation

      IF(K.GE.2.AND.K.LE.4.AND.SWITCH.EQ.0)THEN


!...    Error estimates are to be generated and tested
!...    to see if the tolerance requirements are satisfied.
!
        IFIN = 1
        MSHFLG = 1
        DO  J = 1, MSTAR
          ERREST(J) = 0.D0
        enddo



        YKNEW=1


        DO 60 I = 1, N, 2
         
           X = XI(I)

           CALL APPROX(I,X,Y_MESHNEW(YKNEW),ASAVE(1,5),DUMMY,XI,&
                         N,Z,DMZ,K,NCOMP,MMAX,M,MSTAR,4,DUMMY,0)       


           DO 20 L = 1,MSTAR
             ERR1(L) = DABS(Y_MESHNEW(YKNEW) -  &
             Y_MESHOLD(YKNEW))/2.0d0**(2*k)

             YKNEW = YKNEW + 1

   20      CONTINUE


           IF(I.EQ.(N-1))THEN

              X=XI(I+2)

                CALL APPROX(I+1,X,Y_MESHNEW(YKNEW),ASAVE(1,7),DUMMY,XI,&
                  N,Z,DMZ,K,NCOMP,MMAX,M,MSTAR,4,DUMMY,0)       
      
             
                DO 31 L = 1,MSTAR
                 ERR3(L) =  DABS(Y_MESHNEW(YKNEW) - &
                 Y_MESHOLD(YKNEW))/2.0d0**(2*k)
 
                 YKNEW = YKNEW + 1
   31           CONTINUE
!



           ENDIF


           IF(I.NE.(N-1))THEN

!...       Find component-wise maximum error estimate
!
               DO 40 L = 1,MSTAR
               ERREST(L) = DMAX1(ERREST(L),ERR1(L))
   40          CONTINUE
!

           ELSE



               DO 41 L = 1,MSTAR
               ERREST(L) = DMAX1(ERREST(L),ERR1(L),ERR3(L))
   41          CONTINUE
!

           ENDIF

!...       Test whether the tolerance requirements are satisfied
!...       in the i-th interval.
!
           IF ( IFIN .EQ. 0 )                       GO TO 60
           DO 50 J = 1, NTOL

             LTOLJ = LTOL(J)
             LTJZ = LTOLJ  +  (I-1) * MSTAR
           IF ( ERREST(LTOLJ) .GT. &
                TOLIN(J) * (DABS(Z(LTJZ))+1.D0) )  IFIN = 0

       50  CONTINUE
      60 CONTINUE



            WRITE(IOUT,125)
            LJ = 1

            DO 70 J = 1,NCOMP
              MJ = LJ - 1 + M(J)
              WRITE(IOUT,120) J, (ERREST(L), L= LJ, MJ)
              LJ = MJ + 1
   70       CONTINUE


!     Using SCSA technology for  2 <= k <= 4 immediately after first part of computation with meshpoint error estimation has satisfied the user tolerance (ifin =1)


      IF(IFIN.EQ.1)THEN

!        Set switch = 1 to make sure only doubling is perform to Richardson extrapolation can used going forward.

          SWITCH=1

          IFIN = 0

          MSHFLG=0


         IF(K.EQ.2)THEN
           S_SUPER=4
           S_STAR_SUPER=4
         ENDIF
      
         IF(K.EQ.3)THEN
          S_SUPER=5
          S_STAR_SUPER=6
         ENDIF

        IF(K.EQ.4)THEN
         S_SUPER=6
         S_STAR_SUPER=9
        ENDIF

        ! Set up for SCSA for previous mesh 

        CALL SETUP_PHASE1(NOLD,MESHOLD,ASAVE(1,3),ZOLD,DMZOLD,K,NCOMP,MMAX, &
              M,MSTAR,Y_SUPER,SK_DISCRETE_SUPER,SK_INTERP_SUPER,FSUB)


        KSTORE=1

        DO 156 I = 1, NOLD
         X = MESHOLD(I)

         HD6 = (MESHOLD(I+1) - MESHOLD(I)) / 6.D0
         DO 156 J = 1, 4
           X = X + HD6

           IF ( J.EQ.3 )  X = X + HD6 
    

           CALL EVALUATE(MSTAR,NOLD,MESHOLD,Y_SUPER,X,S_SUPER,S_STAR_SUPER,&
                VALSTR(KSTORE),SK_DISCRETE_SUPER,SK_INTERP_SUPER)

           KSTORE = KSTORE  +  MSTAR
  156   CONTINUE

!
!...  Error estimates are to be generated and tested
!...  to see if the tolerance requirements are satisfied.
!
        IFIN = 1
        MSHFLG = 1

        DO 911 J = 1, MSTAR
  911     ERREST(J) = 0.D0

        CALL SETUP_PHASE1(N,XI,ASAVE(1,3),Z,DMZ,K,NCOMP,MMAX, &
           M,MSTAR,Y_SUPER,SK_DISCRETE_SUPER,SK_INTERP_SUPER,FSUB)


        DO 961 IBACK = 1, N
           I = N + 1 - IBACK
!
!...       The error estimates are obtained by combining values of
!...       the numerical solutions for two meshes.
!...       For each value of Iback=1 to N  we will consider the two
!...       approximations at 2 points in each of
!...       the new subintervals.  we work backwards through
!...       the subinterval so that new values can be stored
!...       in valstr in case they prove to be needed later
!...       for an error estimate. the routine  newmsh
!...       filled in the needed values of the old solution
!...       in valstr.
!
           KNEW = ( 4 * (I-1) + 2 ) * MSTAR + 1
           KSTORE = ( 2 * (I-1) + 1 ) * MSTAR + 1

!          Find a value for z, stored in valstr(knew), at the value x=hi*2/3 on each subinterval
!          Compute the error at this point by finding the weighted difference between valstr(knew) and valstr(kstore)

           X = XI(I) +  (XI(I+1)-XI(I)) * 2.D0 / 3.D0 

           CALL EVALUATE(MSTAR,N,XI,Y_SUPER,X,S_SUPER,S_STAR_SUPER,&
               VALSTR(KNEW),SK_DISCRETE_SUPER,SK_INTERP_SUPER)
             

           DO 921 L = 1,MSTAR
             ERR(L) = WGTERR(L) * DABS(VALSTR(KNEW) -  &
             VALSTR(KSTORE))

             KNEW = KNEW + 1
             KSTORE = KSTORE + 1
   921      CONTINUE

           KNEW = ( 4 * (I-1) + 1 ) * MSTAR + 1
           KSTORE = 2 * (I-1) * MSTAR + 1

!          Find a value for z, stored in valstr(knew), at the value x=xi+hi*1/3 on each subinterval
!          Compute the error at this point by finding the weighted difference between valstr(knew) and valstr(kstore)
!          plus add in the error calculated at the previous point x = xi+hi*2/3

           X = XI(I) +  (XI(I+1)-XI(I)) / 3.D0

           CALL EVALUATE(MSTAR,N,XI,Y_SUPER,X,S_SUPER,S_STAR_SUPER,&
                VALSTR(KNEW),SK_DISCRETE_SUPER,SK_INTERP_SUPER)          


           DO 933 L = 1,MSTAR
             ERR(L) = ERR(L) + WGTERR(L) * DABS(VALSTR(KNEW) - &
             VALSTR(KSTORE))

             KNEW = KNEW + 1
             KSTORE = KSTORE + 1
   933     CONTINUE
!
!...       Find component-wise maximum error estimate
!
           DO 943 L = 1,MSTAR
             ERREST(L) = DMAX1(ERREST(L),ERR(L))
   943      CONTINUE
!

!...       Test whether the tolerance requirements are satisfied
!...       in the i-th interval.
!
           IF ( IFIN .EQ. 0 )                       GO TO 961
           DO 951 J = 1, NTOL
             LTOLJ = LTOL(J)
             LTJZ = LTOLJ  +  (I-1) * MSTAR
           IF ( ERR(LTOLJ) .GT. &
                TOLIN(J) * (DABS(Z(LTJZ))+1.D0) )  IFIN = 0
   951     CONTINUE
   961     CONTINUE


           IF ( IPRINT .GE. 0 )                          RETURN
            WRITE(IOUT,130)
            LJ = 1

            DO 971 J = 1,NCOMP
              MJ = LJ - 1 + M(J)
              WRITE(IOUT,120) J, (ERREST(L), L= LJ, MJ)
              LJ = MJ + 1
   971      CONTINUE

            RETURN

         ENDIF
         ! end of Using SCSA technology for  2 <= k <= 4 immediately after first part of computation has satisfied the user tolerance (ifin =1)

     

!    Using SCSA technology again for  2 <= k <= 4 if first use of SCSA did not satisfy user tolerance. Mesh doubling is still required.


     ELSE IF(K.GE.2.AND.K.LE.4.AND.SWITCH.EQ.1)THEN



       IF(K.EQ.2)THEN
         S_SUPER=4
         S_STAR_SUPER=4
       ENDIF
      
       IF(K.EQ.3)THEN
         S_SUPER=5
         S_STAR_SUPER=6
       ENDIF

       IF(K.EQ.4)THEN
         S_SUPER=6
         S_STAR_SUPER=9
       ENDIF

!
!...  Error estimates are to be generated and tested
!...  to see if the tolerance requirements are satisfied.
!
       IFIN = 1
       MSHFLG = 1
       DO 11 J = 1, MSTAR
   11   ERREST(J) = 0.D0

        CALL SETUP_PHASE1(N,XI,ASAVE(1,3),Z,DMZ,K,NCOMP,MMAX, &
            M,MSTAR,Y_SUPER,SK_DISCRETE_SUPER,SK_INTERP_SUPER,FSUB)


        DO 61 IBACK = 1, N
           I = N + 1 - IBACK
!
!...       The error estimates are obtained by combining values of
!...       the numerical solutions for two meshes.
!...       For each value of Iback=1 to N  we will consider the two
!...       approximations at 2 points in each of
!...       the new subintervals.  we work backwards through
!...       the subinterval so that new values can be stored
!...       in valstr in case they prove to be needed later
!...       for an error estimate. the routine  newmsh
!...       filled in the needed values of the old solution
!...       in valstr.
!
           KNEW = ( 4 * (I-1) + 2 ) * MSTAR + 1
           KSTORE = ( 2 * (I-1) + 1 ) * MSTAR + 1

!          Find a value for z, stored in valstr(knew), at the value x=hi*2/3 on each subinterval
!          Compute the error at this point by finding the weighted difference between valstr(knew) and valstr(kstore)

           X = XI(I) +  (XI(I+1)-XI(I)) * 2.D0 / 3.D0 


           CALL EVALUATE(MSTAR,N,XI,Y_SUPER,X,S_SUPER,S_STAR_SUPER,&
               VALSTR(KNEW),SK_DISCRETE_SUPER,SK_INTERP_SUPER)          

           DO 21 L = 1,MSTAR
             ERR(L) = WGTERR(L) * DABS(VALSTR(KNEW) -  &
             VALSTR(KSTORE))

             KNEW = KNEW + 1
             KSTORE = KSTORE + 1
   21      CONTINUE

           KNEW = ( 4 * (I-1) + 1 ) * MSTAR + 1
           KSTORE = 2 * (I-1) * MSTAR + 1

!          Find a value for z, stored in valstr(knew), at the value x=xi+hi*1/3 on each subinterval
!          Compute the error at this point by finding the weighted difference between valstr(knew) and valstr(kstore)
!          plus add in the error calculated at the previous point x = xi+hi*2/3

           X = XI(I) +  (XI(I+1)-XI(I)) / 3.D0

           CALL EVALUATE(MSTAR,N,XI,Y_SUPER,X,S_SUPER,S_STAR_SUPER,&
                VALSTR(KNEW),SK_DISCRETE_SUPER,SK_INTERP_SUPER)


           DO 33 L = 1,MSTAR
             ERR(L) = ERR(L) + WGTERR(L) * DABS(VALSTR(KNEW) - &
             VALSTR(KSTORE))

             KNEW = KNEW + 1
             KSTORE = KSTORE + 1
   33      CONTINUE
!
!...       Find component-wise maximum error estimate
!
           DO 43 L = 1,MSTAR
             ERREST(L) = DMAX1(ERREST(L),ERR(L))
   43      CONTINUE
!

!...       Test whether the tolerance requirements are satisfied
!...       in the i-th interval.
!
           IF ( IFIN .EQ. 0 )                       GO TO 61
           DO 51 J = 1, NTOL
             LTOLJ = LTOL(J)
             LTJZ = LTOLJ  +  (I-1) * MSTAR
           IF ( ERR(LTOLJ) .GT. &
                TOLIN(J) * (DABS(Z(LTJZ))+1.D0) )  IFIN = 0
   51      CONTINUE
   61      CONTINUE


           IF ( IPRINT .GE. 0 )                          RETURN
           WRITE(IOUT,130)
           LJ = 1
            DO 71 J = 1,NCOMP
             MJ = LJ - 1 + M(J)
             WRITE(IOUT,120) J, (ERREST(L), L= LJ, MJ)
             LJ = MJ + 1
   71      CONTINUE
       RETURN

! Use original colnew code with Richardson extrapolation to test to see if numerical solution satisfies user tolerance


     ELSE ! K=1,5,6,7



      IFIN = 1
      MSHFLG = 1


      DO 19 J = 1, MSTAR
   19   ERREST(J) = 0.D0
      DO 69 IBACK = 1, N
           I = N + 1 - IBACK
!
!...       the error estimates are obtained by combining values of
!...       the numerical solutions for two meshes.
!...       for each value of iback we will consider the two
!...       approximations at 2 points in each of
!...       the new subintervals.  we work backwards through
!...       the subinterval so that new values can be stored
!...       in valstr in case they prove to be needed later
!...       for an error estimate. the routine  newmsh
!...       filled in the needed values of the old solution
!...       in valstr.
!
           KNEW = ( 4 * (I-1) + 2 ) * MSTAR + 1
           KSTORE = ( 2 * (I-1) + 1 ) * MSTAR + 1
           X = XI(I) +  (XI(I+1)-XI(I)) * 2.D0 / 3.D0

           CALL APPROX (I, X, VALSTR(KNEW), ASAVE(1,3), DUMMY, XI,&
                 N, Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 4, DUMMY, 0)

           DO 29 L = 1,MSTAR
             ERR(L) = WGTERR(L) * DABS(VALSTR(KNEW) -   &
            VALSTR(KSTORE))
             KNEW = KNEW + 1
             KSTORE = KSTORE + 1
   29      CONTINUE
           KNEW = ( 4 * (I-1) + 1 ) * MSTAR + 1
           KSTORE = 2 * (I-1) * MSTAR + 1
           X = XI(I) +  (XI(I+1)-XI(I)) / 3.D0

           CALL APPROX (I, X, VALSTR(KNEW), ASAVE(1,2), DUMMY, XI, &
                 N, Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 4, DUMMY, 0)

           DO 39 L = 1,MSTAR
             ERR(L) = ERR(L) + WGTERR(L) * DABS(VALSTR(KNEW) -  &
             VALSTR(KSTORE))
             KNEW = KNEW + 1
             KSTORE = KSTORE + 1
   39      CONTINUE
!
!...       find component-wise maximum error estimate
!
           DO 49 L = 1,MSTAR
             ERREST(L) = DMAX1(ERREST(L),ERR(L))
   49      CONTINUE
!

!...       test whether the tolerance requirements are satisfied
!...       in the i-th interval.
!
           IF ( IFIN .EQ. 0 )                       GO TO 69
           DO 59 J = 1, NTOL
             LTOLJ = LTOL(J)
             LTJZ = LTOLJ  +  (I-1) * MSTAR
           IF ( ERR(LTOLJ) .GT.    &
                TOLIN(J) * (DABS(Z(LTJZ))+1.D0) )  IFIN = 0
   59      CONTINUE
   69      CONTINUE

           IF ( IPRINT .GE. 0 )                          RETURN
            WRITE(IOUT,130)
            LJ = 1
            DO 79 J = 1,NCOMP
             MJ = LJ - 1 + M(J)
             WRITE(IOUT,120) J, (ERREST(L), L= LJ, MJ)
             LJ = MJ + 1
   79      CONTINUE

        RETURN

    ENDIF 
    ! end of first if statement in ERRCHK routine - the decision of what choice of error estimation

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  120 FORMAT (3H U(, I2, 3H) -,4D12.4)
  125 FORMAT (/36H THE MESHPOINT ESTIMATED ERRORS ARE,)
  130 FORMAT (/46H THE CONTINUOUS SOLUTION ESTIMATED ERRORS ARE,)
    

      END SUBROUTINE ERRCHK




!---------------------------------------------------------------------
!                            p a r t  3
!          collocation system setup routines
!---------------------------------------------------------------------
!
      SUBROUTINE LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, &
                 G, W, V, RHS, DMZO, INTEGS, IPVTG, IPVTW, RNORM,&
                 MODE, FSUB, DFSUB, GSUB, DGSUB, GUESS )
!*********************************************************************
!
!   purpose
!         this routine controls the set up and solution of a linear
!      system of collocation equations.
!         the matrix  g  is cast into an almost block diagonal
!      form by an appropriate ordering of the columns and solved
!      using the package of de boor-weiss [5]. the matrix is composed
!      of n blocks. the i-th block has the size
!                  integs(1,i) * integs(2,i).
!      it contains in its last rows the linearized collocation
!      equations, condensed as described in [2],
!      and the linearized side conditions corresponding to
!      the i-th subinterval.  integs(3,i)  steps of gaussian
!      elimination are applied to it to achieve a  partial plu
!      decomposition.  the right hand side vector is put into  rhs
!      and the solution vector is returned in  delz and deldmz.
!
!
!         lsyslv operates according to one of 5 modes:
!      mode = 0 - set up the collocation matrices  v , w , g
!                 and the right hand side  rhs ,  and solve.
!                 (for linear problems only.)
!      mode = 1 - set up the collocation matrices  v , w , g
!                 and the right hand sides  rhs  and  dmzo ,
!                 and solve. also set up  integs .
!                 (first iteration of nonlinear problems only).
!      mode = 2 - set up  rhs  only and compute its norm.
!      mode = 3 - set up  v, w, g  only and solve system.
!      mode = 4 - perform forward and backward substitution only
!                 (do not set up the matrices nor form the rhs).
!
!   variables
!
!      ig,izeta  - pointers to g,zeta respectively
!                       (necessary to keep track of blocks of g
!                       during matrix manipulations)
!      idmz,irhs,iv,iw - pointers to  rhs,v,w rspectively
!      df    - partial derivatives of f from dfsub
!      rnorm - euclidean norm of rhs
!      lside - number of side conditions in current and previous blocks
!      iguess = 1 when current soln is user specified via  guess
!             = 0 otherwise
!
!*********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION  Z(1), DMZ(1), DELZ(1), DELDMZ(1), XI(1), XIOLD(1)
      DIMENSION  G(1), W(1), V(1),  RHS(1) , DMZO(1), DUMMY(1)
      DIMENSION  INTEGS(3,1), IPVTG(1), IPVTW(1)
      DIMENSION  ZVAL(40), F(40), DGZ(40), DMVAL(20), DF(800), AT(28)
!

      COMMON /COLOUT/ PRECIS, IOUT, IPRINT
      COMMON /COLLOC/ RHO(7), COEF(49)
      COMMON /COLORD/ K, NCOMP, MSTAR, KD,  MMAX, M(20)
      COMMON /COLSID/ ZETA(40), ALEFT, ARIGHT, IZETA, IZSAVE
      COMMON /COLAPR/ N, NOLD, NMAX, NZ, NDMZ
      COMMON /COLNLN/ NONLIN, ITER, LIMIT, ICARE, IGUESS
      COMMON /COLBAS/ B(28), ACOL(28,7), ASAVE(28,7)
!
      EXTERNAL DFSUB, DGSUB, FSUB, GSUB, GUESS
      OPTIONAL :: DFSUB, DGSUB

!-----------------------------------------------------------------------------------------
!
!     Common Blocks
!
!     colout- 1) Precis is a constant = 100*macine unit roundoff
!             2) Iout is a constant = 6 for writing output to the screen
!             3) Iprint is equal to ipar(7) and therefore iprint
!         
!                = -1 for full diagnostic printout
!                = 0 for selected printout
!                = 1 for no printout
!      
!     colloc - 1) Rho are the k collocation points on (0,1)
!              2) Coef are the precomputed derivatives of the basis
!
!     colord - 1) K is the number of collocation points
!              2) Nc is equal to ncomp is the number of differential equations
!              3) Mstar is the sum of the orders of the differential equations
!                 which is equal to the number of solution and derivative tolerances
!                 which is also equal to the number of differential equations in an
!                 equivalent first order system.
!              4) Kd is the number of collocation points multiplied by the number of
!                 differential equations.
!              5) Mmax is the maxmium order of the set of differential equations
!              6) M is an array which stores all of the orders of each differential equation
!
!    colapr - 1) N is the number of subintervals in the current mesh
!             2) Nold is the number of subintervals in the previous mesh
!             3) Nmax is the maximum number of subintervals allowed
!             4) Nz is the number of components in the Z solution vector equal to
!                Mstar*(N+1)
!             5) Ndmz is equal to the number of collocation points multiplied by the number of differential
!                equations multiplied by the number of subintervals.
!                This is equal to the number of components of Dmz equal to Kd*N
!
!   colsid - 1) tzeta is an array equal to zeta(i) so
!
!               tzeta(j) - j-th side condition point (boundary point). must
!                          have  tzeta(j) .le. tzeta(j+1). all side condition
!                          points must be mesh points in all meshes used,
!                          see description of ipar(11) and fixpnt below. 
!
!               (this variable name is not used in computation)
!
!             2) tleft is equal to aleft the left boundary point
!                (this variable name is not used in computation)
!             3) tright is equal to aright the right boundary point
!                (this variable name is not used in computation)
!             4) izeta points to zeta (allows incrementation by 1)
!             5) idum equals izave equals izeta
!
!  colnln    1) nonlin is equal to ipar(1)
!                = 0 if the problem is linear
!                = 1 if the problem is nonlinear
!
!            2) iter is the number of iterations (counted only when Jacobian inversions are performed)
!            3) limit is the maximum number of newton iterations per mesh
!            4) icare = ipar(10) initially
!
!               = 0 if the problem is regular
!               = 1 if the first relax factor is =rstart, and the
!                 nonlinear iteration does not rely on past covergence
!                 (use for an extra sensitive nonlinear problem only).
!               = 2 if we are to return immediately upon  (a) two
!                 successive nonconvergences, or  (b) after obtaining
!                 error estimate for the first time.
!
!               Later in the program the meaning changes slightly to
!
!               = -1 no convergence occured for a reular problem
!               = 0 a regular problem
!               = 1 a sensitive problem
!               = 2 used for continuation ( see description of ipar(10)
!                 in colnew).
!             5) iguess = ipar(9) 
!                 = 0 if no initial guess for the solution is
!                 provided.
!                 = 1 if an initial guess is provided by the user
!                 in subroutine  guess.
!                 = 2 if an initial mesh and approximate solution
!                 coefficients are provided by the user in  fspace.
!                 (the former and new mesh are the same).
!                 = 3 if a former mesh and approximate solution
!                 coefficients are provided by the user in fspace,
!                 and the new mesh is to be taken twice as coarse;
!                 i.e.,every second point from the former mesh.
!                 = 4 if in addition to a former initial mesh and
!                 approximate solution coefficients, a new mesh
!                 is provided in fspace as well.
!                 (see description of output for further details
!                 on iguess = 2, 3, and 4. see also simple continuation documentation above)
!
!       colbas - 1) B=
!                   = rkb - the runge kutta basis (0th to (m-1)th derivatives)
!                   computed in rkbas at S (the relative position for which the basis is to be evaluated)
!                   (0.le.S.le.1)
!                     
!                2) acol is the runge kutta coefficients at the collocation points (rho) = S
!                3) asave is the runge-kutta coefficients at (1/6,1/3,2/3,5/6) =S used in
!                   newmesh and errchk 
!   
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Contrl
!               
!     Calls to : Approx, Sbblok, Dmzsol, Gderiv, Vwblok, Gblock, Fcblok
!
!-----------------------------------------------------------------------------------------------------------------------

!     For linear case - mode = 0 so M1 = 1, therefore GOTO 10      

      M1 = MODE + 1
      GO TO (10, 30, 30, 30, 310), M1
!
!...  Linear problem initialization
!
   10 DO 20 I=1,MSTAR
   20 ZVAL(I) = 0.D0
!
!...  Initialization - 
!
!      ig,izeta  - pointers to g,zeta respectively
!                       (necessary to keep track of blocks of g
!                       during matrix manipulations)
!      idmz,irhs,iv,iw - pointers to idmz,rhs,v,w rspectively

   30 IDMZ = 1
      IDMZO = 1
      IRHS = 1
      IG = 1
      IW = 1
      IV = 1
      IZETA = 1
      LSIDE = 0
      IOLD = 1
      NCOL = 2 * MSTAR
      RNORM = 0.D0

!     If mode is greater than one no need to set up integs. GOTO 80

      IF ( MODE .GT. 1 )                            GO TO 80
!
!...  Build integs (describing block structure of matrix)
!
      DO 70 I = 1,N
           INTEGS(2,I) = NCOL
           IF (I .LT. N)                            GO TO 40
           INTEGS(3,N) = NCOL
           LSIDE = MSTAR
           GO TO 60
   40      INTEGS(3,I) = MSTAR
   50      IF( LSIDE .EQ. MSTAR )                   GO TO 60
           IF ( ZETA(LSIDE+1) .GE. XI(I)+PRECIS )   GO TO 60
           LSIDE = LSIDE + 1
           GO TO 50
   60      NROW = MSTAR + LSIDE
   70 INTEGS(1,I) = NROW
   80 CONTINUE

!     If setting up rhs only and computing its norm GOTO 90

      IF ( MODE .EQ. 2 )                            GO TO 90
!
!...  Zero the matrices to be computed
!
      LW = KD * KD * N
      DO 84 L = 1, LW
   84   W(L) = 0.D0
!
!...  The do loop 290 sets up the linear system of equations.
!
  90  CONTINUE


      DO 290 I=1, N
!
!...       Construct a block of  A  and a corresponding piece of  rhs.
!
           XII = XI(I)
           H = XI(I+1) - XI(I)
           NROW = INTEGS(1,I)
!
!...       Go thru the ncomp collocation equations and side conditions
!...       in the i-th subinterval
!
  100      IF ( IZETA .GT. MSTAR )                  GO TO 140
           IF ( ZETA(IZETA) .GT. XII + PRECIS )      GO TO 140
!
!...       Build equation for a side condition.
!
!          If linear problem

           IF ( MODE .EQ. 0 )                       GO TO 110
           IF ( IGUESS .NE. 1 )                     GO TO 102
!
!...       Case where user provided current approximation


           CALL GUESS (XII, ZVAL, DMVAL)


           GO TO 110
!
!...       other nonlinear case
!
  102      IF ( MODE .NE. 1 )                       GO TO 106


           CALL APPROX (IOLD, XII, ZVAL, AT, COEF, XIOLD, NOLD, Z, DMZ, &
                        K, NCOMP, MMAX, M, MSTAR, 2, DUMMY, 0)

!          ENDIF

           GO TO 110


  106      CALL APPROX (I, XII, ZVAL, AT, DUMMY, XI, N, Z, DMZ, &
                        K, NCOMP, MMAX, M, MSTAR, 1, DUMMY, 0)
  108      IF ( MODE .EQ. 3 )                       GO TO 120
!
!...       find  rhs  boundary value.
!


  110      CALL GSUB (IZETA, ZVAL, GVAL)


           RHS(NDMZ+IZETA) = -GVAL
           RNORM = RNORM + GVAL**2
           IF ( MODE .EQ. 2 )                       GO TO 130
!
!...       build a row of  a  corresponding to a boundary point
!


  120      CALL GDERIV (G(IG), NROW, IZETA, ZVAL, DGZ, 1, GSUB, DGSUB)
  130      IZETA = IZETA + 1
           GO TO 100
!
!...       assemble collocation equations
!
  140      DO 220 J = 1, K
             HRHO = H * RHO(J)
             XCOL = XII + HRHO
!
!...         this value corresponds to a collocation (interior)
!...         point. build the corresponding  ncomp  equations.
!
             IF ( MODE .EQ. 0 )                     GO TO 200
             IF ( IGUESS .NE. 1 )                   GO TO 160
!
!...         use initial approximation provided by the user.
!
             CALL GUESS (XCOL, ZVAL, DMZO(IRHS) )
             GO TO 170
!
!...         find  rhs  values
!
  160        IF ( MODE .NE. 1 )                     GO TO 190


             CALL APPROX (IOLD, XCOL, ZVAL, AT, COEF, XIOLD, NOLD,&
                  Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 2, DMZO(IRHS), 1)
!
  170        CALL FSUB (XCOL, ZVAL, F)
             DO 180 JJ = 1, NCOMP
               VALUE = DMZO(IRHS) - F(JJ)
               RHS(IRHS) = - VALUE
               RNORM = RNORM + VALUE**2
               IRHS = IRHS + 1
  180        CONTINUE
             GO TO 210
!
!...         evaluate former collocation solution
!


  190        CALL APPROX (I, XCOL, ZVAL, ACOL(1,J), COEF, XI, N, &
                  Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 4, DUMMY, 0)
             IF ( MODE .EQ. 3 )                     GO TO 210
!
!...         fill in  rhs  values (and accumulate its norm).
!
             CALL FSUB (XCOL, ZVAL, F)
             DO 195 JJ = 1, NCOMP
               VALUE = DMZ(IRHS) - F(JJ)
               RHS(IRHS) = - VALUE
               RNORM = RNORM + VALUE**2
               IRHS = IRHS + 1
  195        CONTINUE
             GO TO 220
!
!...         the linear case
!
  200        CALL FSUB (XCOL, ZVAL, RHS(IRHS))
             IRHS = IRHS + NCOMP
!
!...         fill in ncomp rows of  w and v
!


  210        CALL VWBLOK (XCOL, HRHO, J, W(IW), V(IV), IPVTW(IDMZ), &
             KD, ZVAL, DF, ACOL(1,J), DMZO(IDMZO), NCOMP, FSUB, DFSUB, MSING)


             IF ( MSING .NE. 0 )                    RETURN
  220      CONTINUE
!
!...       build global bvp matrix  g
!
           IF ( MODE .NE. 2 ) &
            CALL GBLOCK (H, G(IG), NROW, IZETA, W(IW), V(IV), KD, &
                        DUMMY, DELDMZ(IDMZ), IPVTW(IDMZ), 1 )
           IF ( I .LT. N )                          GO TO 280
           IZSAVE = IZETA
  240      IF ( IZETA .GT. MSTAR )                  GO TO 290
!
!...       build equation for a side condition.
!
           IF ( MODE .EQ. 0 )                       GO TO 250
           IF ( IGUESS .NE. 1 )                     GO TO 245
!
!...       case where user provided current approximation
!
           CALL GUESS (ARIGHT, ZVAL, DMVAL)
           GO TO 250
!
!...       other nonlinear case
!
  245      IF ( MODE .NE. 1 )                       GO TO 246



           CALL APPROX (NOLD+1, ARIGHT, ZVAL, AT, COEF, XIOLD, NOLD, &
                Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 1, DUMMY, 0)
           GO TO 250


  246      CALL APPROX (N+1, ARIGHT, ZVAL, AT, COEF, XI, N, &
                Z, DMZ, K, NCOMP, MMAX, M, MSTAR, 1, DUMMY, 0)
  248      IF ( MODE .EQ. 3 )                       GO TO 260
!
!...       find  rhs  boundary value.
!
  250      CALL GSUB (IZETA, ZVAL, GVAL)
           RHS(NDMZ+IZETA) = - GVAL
           RNORM = RNORM + GVAL**2
           IF ( MODE .EQ. 2 )                       GO TO 270
!
!...       build a row of  a  corresponding to a boundary point
!
  260      CALL GDERIV (G(IG), NROW, IZETA+MSTAR, ZVAL, DGZ, 2, GSUB, DGSUB)
  270      IZETA = IZETA + 1
           GO TO 240
!
!...       update counters -- i-th block completed
!
  280      IG = IG + NROW * NCOL
           IV = IV + KD * MSTAR
           IW = IW + KD * KD
           IDMZ = IDMZ + KD
           IF ( MODE .EQ. 1 )  IDMZO = IDMZO + KD
  290 CONTINUE
!
!...       assembly process completed
!
      IF ( MODE .EQ. 0 .OR. MODE .EQ. 3 )           GO TO 300
      RNORM = DSQRT(RNORM / DFLOAT(NZ+NDMZ))
      IF ( MODE .NE. 2 )                            GO TO 300
      RETURN
!
!...  solve the linear system.
!
!...  matrix decomposition
!
  300 CALL FCBLOK (G, INTEGS, N, IPVTG, DF, MSING)
!
!...  check for singular matrix
!
      MSING = - MSING
      IF( MSING .NE. 0 )                            RETURN
!
!...  perform forward and backward substitution for mode=4 only.
!
  310 CONTINUE
      DO 311 L = 1, NDMZ
        DELDMZ(L) = RHS(L)
  311 CONTINUE
      IZ = 1
      IDMZ = 1
      IW = 1
      IZET = 1
      DO 320 I=1, N
         NROW = INTEGS(1,I)
         IZETA = NROW + 1 - MSTAR
         IF ( I .EQ. N ) IZETA = IZSAVE
  322    IF ( IZET .EQ. IZETA )                     GO TO 324
           DELZ(IZ-1+IZET) = RHS(NDMZ+IZET)
           IZET = IZET + 1
         GO TO 322
  324    H = XI(I+1) - XI(I)
         CALL GBLOCK (H, G(1), NROW, IZETA, W(IW), V(1), KD, &
                      DELZ(IZ), DELDMZ(IDMZ), IPVTW(IDMZ), 2 )
         IZ = IZ + MSTAR
         IDMZ = IDMZ + KD
         IW = IW + KD * KD
         IF ( I .LT. N )                            GO TO 320
  326    IF ( IZET .GT. MSTAR )                     GO TO 320
           DELZ(IZ-1+IZET) = RHS(NDMZ+IZET)
           IZET = IZET + 1
         GO TO 326
  320 CONTINUE
!
!...  perform forward and backward substitution for mode=0,2, or 3.
!
      CALL SBBLOK (G, INTEGS, N, IPVTG, DELZ)
!
!...  finaly find deldmz
!
      CALL DMZSOL (KD, MSTAR, N, V, DELZ, DELDMZ)
!
      IF ( MODE .NE. 1 )                            RETURN
      DO 321 L = 1, NDMZ
        DMZ(L) = DMZO(L)
  321 CONTINUE
      IZ = 1
      IDMZ = 1
      IW = 1
      IZET = 1
      DO 350 I=1, N
         NROW = INTEGS(1,I)
         IZETA = NROW + 1 - MSTAR
         IF ( I .EQ. N ) IZETA = IZSAVE
  330    IF ( IZET .EQ. IZETA )                     GO TO 340
           Z(IZ-1+IZET) = DGZ(IZET)
           IZET = IZET + 1
         GO TO 330
  340    H = XI(I+1) - XI(I)
         CALL GBLOCK (H, G(1), NROW, IZETA, W(IW), DF, KD, &
                      Z(IZ), DMZ(IDMZ), IPVTW(IDMZ), 2 )
         IZ = IZ + MSTAR
         IDMZ = IDMZ + KD
         IW = IW + KD * KD
         IF ( I .LT. N )                            GO TO 350
  342    IF ( IZET .GT. MSTAR )                     GO TO 350
            Z(IZ-1+IZET) = DGZ(IZET)
            IZET = IZET + 1
         GO TO 342
  350 CONTINUE
      CALL SBBLOK (G, INTEGS, N, IPVTG, Z)
!
!...  finaly find dmz
!
      CALL DMZSOL (KD, MSTAR, N, V, Z, DMZ)
!
      RETURN


      END SUBROUTINE LSYSLV



      SUBROUTINE GDERIV ( GI, NROW, IROW, ZVAL, DGZ, MODE, GSUB, DGSUB)
!
!**********************************************************************
!
!   purpose:
!
!      construct a collocation matrix row according to mode:
!      mode = 1  -  a row corresponding to a initial condition
!                   (i.e. at the left end of the subinterval).
!      mode = 2  -  a row corresponding to a final condition.
!
!   variables:
!
!      gi     - the sub-block of the global bvp matrix in
!               which the equations are to be formed.
!      nrow   - no. of rows in gi.
!      irow   - the row in gi to be used for equations.
!      zval   - z(xi)
!      dg     - the derivatives of the side condition.
!
!**********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION GI(NROW,1), ZVAL(1), DGZ(1), DG(40)
!
      COMMON /COLORD/ KDUM, NDUM, MSTAR, KD, MMAX, M(20)
      COMMON /COLSID/ ZETA(40), ALEFT, ARIGHT, IZETA, IDUM
      COMMON /COLNLN/ NONLIN, ITER, LIMIT, ICARE, IGUESS

      EXTERNAL DGSUB,GSUB
      OPTIONAL :: DGSUB

!---------------------------------------------------------------------------------
!
!     Common Blocks
!
!     colord - 1) K is number of collocation points
!
!              2) Nc is equal to ncomp is the number of differential equations
!              3) Mstar is the sum of the orders of the differential equations
!                 which is equal to the number of solution and derivative tolerances
!                 which is also equal to the number of differential equations in an
!                 equivalent first order system.
!              4) Kd is the number of collocation points multiplied by the number of
!                 differential equations.
!              5) Mmax is the maxmium order of the set of differential equations
!              6) M is an array which stores all of the orders of each differential equation
!
!   colsid - 1) tzeta is an array equal to zeta(i) so
!
!               tzeta(j) - j-th side condition point (boundary point). must
!                          have  tzeta(j) .le. tzeta(j+1). all side condition
!                          points must be mesh points in all meshes used,
!                          see description of ipar(11) and fixpnt below. 
!
!               (this variable name is not used in computation)
!
!             2) tleft is equal to aleft the left boundary point
!                (this variable name is not used in computation)
!             3) tright is equal to aright the right boundary point
!                (this variable name is not used in computation)
!             4) izeta points to zeta (allows incrementation by 1)
!             5) idum equals izave equals izeta
!
!  colnln    1) nonlin is equal to ipar(1)
!                = 0 if the problem is linear
!                = 1 if the problem is nonlinear
!
!            2) iter is the number of iterations (counted only when Jacobian inversions are performed)
!            3) limit is the maximum number of newton iterations per mesh
!            4) icare = ipar(10) initially
!
!               = 0 if the problem is regular
!               = 1 if the first relax factor is =rstart, and the
!                 nonlinear iteration does not rely on past covergence
!                 (use for an extra sensitive nonlinear problem only).
!               = 2 if we are to return immediately upon  (a) two
!                 successive nonconvergences, or  (b) after obtaining
!                 error estimate for the first time.
!
!               Later in the program the meaning changes slightly to
!
!               = -1 no convergence occured for a reular problem
!               = 0 a regular problem
!               = 1 a sensitive problem
!               = 2 used for continuation ( see description of ipar(10)
!                 in colnew).
!             5) iguess = ipar(9) 
!                 = 0 if no initial guess for the solution is
!                 provided.
!                 = 1 if an initial guess is provided by the user
!                 in subroutine  guess.
!                 = 2 if an initial mesh and approximate solution
!                 coefficients are provided by the user in  fspace.
!                 (the former and new mesh are the same).
!                 = 3 if a former mesh and approximate solution
!                 coefficients are provided by the user in fspace,
!                 and the new mesh is to be taken twice as coarse;
!                 i.e.,every second point from the former mesh.
!                 = 4 if in addition to a former initial mesh and
!                 approximate solution coefficients, a new mesh
!                 is provided in fspace as well.
!                 (see description of output for further details
!                 on iguess = 2, 3, and 4. see also simple continuation documentation above)
!
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Lsyslv
!               
!     Calls to : No subroutines
!
!-----------------------------------------------------------------------------------------------------------------------
!
!...  zero jacobian dg
!
      DO 10 J=1,MSTAR
   10   DG(J) = 0.D0
!
!...  evaluate jacobian dg
!

      IF(PRESENT(DGSUB)) THEN
       
        CALL DGSUB (IZETA, ZVAL, DG)

      ELSE

        CALL DBCDY_DIVDIFF(IZETA,ZVAL,MSTAR,GSUB,DG)

      ENDIF

!
!...  evaluate  dgz = dg * zval  once for a new mesh
!
      IF (NONLIN .EQ. 0 .OR. ITER .GT. 0)           GO TO 30
      DOT = 0.D0
      DO 20 J = 1, MSTAR
   20   DOT = DOT  +  DG(J) * ZVAL(J)
      DGZ(IZETA) = DOT
!
!...  branch according to  m o d e
!
   30 IF ( MODE .EQ. 2 )                            GO TO 50
!
!...  provide coefficients of the j-th linearized side condition.
!...  specifically, at x=zeta(j) the j-th side condition reads
!...  dg(1)*z(1) + ... +dg(mstar)*z(mstar) + g = 0
!
!
!...  handle an initial condition
!
      DO 40 J = 1, MSTAR
        GI(IROW,J) =  DG(J)
   40 GI(IROW,MSTAR+J) = 0.D0
      RETURN
!
!...  handle a final condition
!
   50 DO 60 J= 1, MSTAR
        GI(IROW,J) = 0.D0
   60 GI(IROW,MSTAR+J) = DG(J)
      RETURN


      END SUBROUTINE GDERIV



      SUBROUTINE VWBLOK (XCOL, HRHO, JJ, WI, VI, IPVTW, KD, ZVAL,&
                         DF, ACOL, DMZO, NCOMP, FSUB, DFSUB, MSING)
!
!**********************************************************************
!
!   purpose:
!
!      construct a group of  ncomp  rows of the matrices  wi  and  vi.
!      corresponding to an interior collocation point.
!
!
!   variables:
!
!      xcol   - the location of the collocation point.
!      jj     - xcol is the jj-th of k collocation points
!               in the i-th subinterval.
!      wi,vi  - the i-th block of the collocation matrix
!               before parameter condensation.
!      kd     - no. of rows in vi and wi .
!      zval   - z(xcol)
!      df     - the jacobian at xcol .
!      jcomp  - counter for the component being dealt with.
!
!**********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION WI(KD,1), VI(KD,1), ZVAL(1), DMZO(1), DF(NCOMP,1)
      DIMENSION IPVTW(1),  HA(7,4), ACOL(7,4), BASM(5)
!
      COMMON /COLORD/ K, NCDUM, MSTAR, KDUM, MMAX, M(20)
      COMMON /COLNLN/ NONLIN, ITER, LIMIT, ICARE, IGUESS

      EXTERNAL FSUB, DFSUB
      OPTIONAL :: DFSUB

!------------------------------------------------------------------------
!
!     Common Blocks
!
!
!     colord - 1) K is the number of collocation points
!              2) Nc is equal to ncomp is the number of differential equations
!              3) Mstar is the sum of the orders of the differential equations
!                 which is equal to the number of solution and derivative tolerances
!                 which is also equal to the number of differential equations in an
!                 equivalent first order system.
!              4) Kd is the number of collocation points multiplied by the number of
!                 differential equations.
!              5) Mmax is the maxmium order of the set of differential equations
!              6) M is an array which stores all of the orders of each differential equation
!
!
!  colnln    1) nonlin is equal to ipar(1)
!                = 0 if the problem is linear
!                = 1 if the problem is nonlinear
!
!            2) iter is the number of iterations (counted only when Jacobian inversions are performed)
!            3) limit is the maximum number of newton iterations per mesh
!            4) icare = ipar(10) initially
!
!               = 0 if the problem is regular
!               = 1 if the first relax factor is =rstart, and the
!                 nonlinear iteration does not rely on past covergence
!                 (use for an extra sensitive nonlinear problem only).
!               = 2 if we are to return immediately upon  (a) two
!                 successive nonconvergences, or  (b) after obtaining
!                 error estimate for the first time.
!
!               Later in the program the meaning changes slightly to
!
!               = -1 no convergence occured for a reular problem
!               = 0 a regular problem
!               = 1 a sensitive problem
!               = 2 used for continuation ( see description of ipar(10)
!                 in colnew).
!             5) iguess = ipar(9) 
!                 = 0 if no initial guess for the solution is
!                 provided.
!                 = 1 if an initial guess is provided by the user
!                 in subroutine  guess.
!                 = 2 if an initial mesh and approximate solution
!                 coefficients are provided by the user in  fspace.
!                 (the former and new mesh are the same).
!                 = 3 if a former mesh and approximate solution
!                 coefficients are provided by the user in fspace,
!                 and the new mesh is to be taken twice as coarse;
!                 i.e.,every second point from the former mesh.
!                 = 4 if in addition to a former initial mesh and
!                 approximate solution coefficients, a new mesh
!                 is provided in fspace as well.
!                 (see description of output for further details
!                 on iguess = 2, 3, and 4. see also simple continuation documentation above)
!
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Lsyslv
!               
!     Calls to : No subroutines
!
!-----------------------------------------------------------------------------------------------------------------------

!
!...  if jj = 1 initialize  wi .
!
      IF ( JJ .GT. 1 )                              GO TO 30
      DO 10 ID = 1, KD
        WI(ID,ID) = 1.D0
   10 CONTINUE
!
!...  calculate local basis
!
   30        FACT = 1.D0
             DO 150 L=1,MMAX
                FACT = FACT * HRHO / DFLOAT(L)
                BASM(L) = FACT
                DO 150 J=1,K
                   HA(J,L) = FACT * ACOL(J,L)
  150        CONTINUE
!
!... zero jacobian
!
      DO 40 JCOL = 1, MSTAR
        DO 40 IR = 1, NCOMP
   40 DF(IR,JCOL) = 0.D0
!
!...  build ncomp rows for interior collocation point x.
!...  the linear expressions to be constructed are:
!...   (m(id))
!...  u     -  df(id,1)*z(1) - ... - df(id,mstar)*z(mstar)
!...   id
!...  for id = 1 to ncomp.

      IF(PRESENT(DFSUB)) THEN
       
        CALL DFSUB (XCOL, ZVAL, DF)

      ELSE

        CALL DFDY_DIVDIFF(XCOL,ZVAL,NCOMP,MSTAR,FSUB,DF)

      ENDIF

!
      I0 = (JJ-1) * NCOMP
      I1 = I0 + 1
      I2 = I0 + NCOMP
!
!...  evaluate  dmzo = dmz - df * zval  once for a new mesh
!
      IF (NONLIN .EQ. 0 .OR. ITER .GT. 0)          GO TO 60
      DO 50 J = 1, MSTAR
        FACT = - ZVAL(J)
        DO 50 ID = 1, NCOMP
          DMZO(I0+ID) = DMZO(I0+ID)  +  FACT * DF(ID,J)
  50  CONTINUE
!
!...  loop over the  ncomp  expressions to be set up for the
!...  current collocation point.
!
   60 DO 70 J = 1, MSTAR
        DO 70 ID = 1, NCOMP
          VI(I0+ID,J) = DF(ID,J)
   70 CONTINUE
      JN = 1
      DO 140 JCOMP = 1, NCOMP
         MJ = M(JCOMP)
         JN = JN + MJ
         DO 130 L = 1, MJ
            JV = JN - L
            JW = JCOMP
            DO 90 J = 1, K
              AJL = - HA(J,L)
              DO 80 IW = I1, I2
                 WI(IW,JW) = WI(IW,JW)  +  AJL * VI(IW,JV)
   80         CONTINUE
   90       JW = JW + NCOMP
            LP1 = L + 1
            IF ( L .EQ. MJ )                        GO TO 130
            DO 110 LL = LP1, MJ
              JDF = JN - LL
              BL = BASM(LL-L)
              DO 100 IW = I1, I2
                VI(IW,JV) = VI(IW,JV)  +  BL * VI(IW,JDF)
  100         CONTINUE
  110       CONTINUE
  130    CONTINUE
  140 CONTINUE
      IF ( JJ .LT. K )                          RETURN
!
!   ...decompose the wi block and solve for the mstar columns of vi
!
!
!...  do parameter condensation
!
      MSING = 0
      CALL DGEFA  (WI, KD, KD, IPVTW, MSING)
!
!...   check for singularity
!
      IF ( MSING .NE. 0 )                         RETURN
      DO 250 J= 1,MSTAR
         CALL DGESL  (WI, KD, KD, IPVTW, VI(1,J), 0)
  250 CONTINUE
      RETURN


      END SUBROUTINE VWBLOK

      SUBROUTINE GBLOCK (H, GI, NROW, IROW, WI, VI, KD, &
                         RHSZ, RHSDMZ, IPVTW, MODE)
!
!**********************************************************************
!
!   purpose:
!
!      construct collocation matrix rows according to mode:
!      mode = 1  -  a group of  mstar    rows corresponding
!                   an interior mesh interval.
!           = 2  -  corresponding right hand side
!
!   variables:
!
!      h      - the  local stepsize.
!      gi     - the sub-block of the collocation matrix in
!               which the equations are to be formed.
!      wi     - the sub-block of noncondensed collocation equations,
!               left-hand side part.
!      vi     - the sub-block of noncondensed collocation equations,
!               right-hand side part.
!      rhsdmz - the inhomogenous term of the uncondensed collocation
!               equations.
!      rhsz   - the inhomogenous term of the condensed collocation
!               equations.
!      nrow   - no. of rows in gi.
!      irow   - the first row in gi to be used for equations.
!
!**********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION HB(7,4), BASM(5)
      DIMENSION GI(NROW,1), WI(1), VI(KD,1)
      DIMENSION RHSZ(1), RHSDMZ(1), IPVTW(1)
!
      COMMON /COLORD/  K, NCOMP, MSTAR, KDUM, MMAX, M(20)
      COMMON /COLBAS/ B(7,4), ACOL(28,7), ASAVE(28,7)
!
!-------------------------------------------------------------------------------------
!
!     Common Blocks
!
!
!     colord - 1) K is the number of collocation points
!              2) Nc is equal to ncomp is the number of differential equations
!              3) Mstar is the sum of the orders of the differential equations
!                 which is equal to the number of solution and derivative tolerances
!                 which is also equal to the number of differential equations in an
!                 equivalent first order system.
!              4) Kd is the number of collocation points multiplied by the number of
!                 differential equations.
!              5) Mmax is the maxmium order of the set of differential equations
!              6) M is an array which stores all of the orders of each differential equation
!
!       colbas - 1) B=
!                   = rkb - the runge kutta basis (0th to (m-1)th derivatives)
!                   computed in rkbas at S (the relative position for which the basis is to be evaluated)
!                   (0.le.S.le.1)
!                     
!                2) acol is the runge kutta coefficients at the collocation points (rho) = S
!                3) asave is the runge-kutta coefficients at (1/6,1/3,2/3,5/6) =S used in
!                   newmesh and errchk 
!   
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Lsyslv
!               
!     Calls to : No subroutines
!
!-----------------------------------------------------------------------------------------------------------------------

!...  compute local basis
!
      FACT = 1.D0
      BASM(1) = 1.D0
      DO 30 L=1,MMAX
         FACT = FACT * H / DFLOAT(L)
         BASM(L+1) = FACT
         DO 20 J=1,K
   20       HB(J,L) = FACT * B(J,L)
   30 CONTINUE
!
!...  branch according to  m o d e
!
      GO TO (40, 110), MODE
!
!...  set right gi-block to identity
!
   40 CONTINUE
      DO 60 J = 1, MSTAR
        DO 50 IR = 1, MSTAR
          GI(IROW-1+IR,J) = 0.D0
   50   GI(IROW-1+IR,MSTAR+J) = 0.D0
   60 GI(IROW-1+J,MSTAR+J) = 1.D0
!
!...  compute the block gi
!
      IR = IROW
      DO 100 ICOMP = 1, NCOMP
         MJ = M(ICOMP)
         IR = IR + MJ
         DO 90 L = 1, MJ
            ID = IR - L
            DO 80 JCOL = 1, MSTAR
               IND = ICOMP
               RSUM = 0.D0
               DO 70 J = 1, K
                  RSUM = RSUM  -  HB(J,L) * VI(IND,JCOL)
   70          IND = IND + NCOMP
               GI(ID,JCOL) = RSUM
   80       CONTINUE
            JD = ID - IROW
            DO 85 LL = 1, L
               GI(ID,JD+LL) = GI(ID,JD+LL) - BASM(LL)
   85       CONTINUE
   90    CONTINUE
  100 CONTINUE
      RETURN
!
!...  compute the appropriate piece of  rhsz
!
  110 CONTINUE
      CALL DGESL  (WI, KD, KD, IPVTW, RHSDMZ, 0)
      IR = IROW
      DO 140 JCOMP = 1, NCOMP
         MJ = M(JCOMP)
         IR = IR + MJ
         DO 130 L = 1, MJ
            IND = JCOMP
            RSUM = 0.D0
            DO 120 J = 1, K
               RSUM = RSUM  +  HB(J,L) * RHSDMZ(IND)
  120       IND = IND + NCOMP
            RHSZ(IR-L) = RSUM
  130    CONTINUE
  140 CONTINUE
      RETURN


      END SUBROUTINE GBLOCK



!----------------------------------------------------------------------
!                             p a r t  4
!               polynomial and service routines
!----------------------------------------------------------------------
!
      SUBROUTINE APPSLN (X, Z, FSPACE, ISPACE,FSUB)
!
!*****************************************************************
!
!     purpose
!
!           set up a standard call to  approx  to evaluate the
!           approximate solution  z = z( u(x) ) using SCSA technology at a point x
!           (it has been computed by a call to  colnew ).
!           the parameters needed for  approx  are retrieved
!           from the work arrays  ispace  and  fspace .
!
!*****************************************************************
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Main Program
!               
!     Calls to : Approx
!
!-----------------------------------------------------------------------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION Z(1), FSPACE(*), ISPACE(*), A(28),DUMMY(1)
     
      
      DOUBLE PRECISION Y_SUPER(160*(3000+1))
      DOUBLE PRECISION SK_DISCRETE_SUPER(3000*10*160)
      DOUBLE PRECISION SK_INTERP_SUPER(3000*10*160)

      INTEGER S_SUPER
      INTEGER S_STAR_SUPER
!
      INTEGER IMAXE

      COMMON /APPCOL/IMAXE  

      EXTERNAL FSUB

      IMAXE=314

      IF(ISPACE(2).EQ.2)THEN
        S_SUPER=4
        S_STAR_SUPER=4
      ENDIF
      
      IF(ISPACE(2).EQ.3)THEN
        S_SUPER=5
        S_STAR_SUPER=6
      ENDIF

      IF(ISPACE(2).EQ.4)THEN
        S_SUPER=6
        S_STAR_SUPER=9
      ENDIF

      IS6 = ISPACE(6)
      IS5 = ISPACE(1) + 2
      IS4 = IS5 + ISPACE(4) * (ISPACE(1) + 1)
      I=1


      IF(ISPACE(2).GE.2.AND.ISPACE(2).LE.4)THEN

        CALL SETUP_PHASE1(ISPACE(1),FSPACE(1),A,FSPACE(IS5),FSPACE(IS4), &
           ISPACE(2),ISPACE(3),ISPACE(5),ISPACE(8),ISPACE(4), &
           Y_SUPER,SK_DISCRETE_SUPER,SK_INTERP_SUPER,FSUB)
     
  
        CALL EVALUATE(ISPACE(4),ISPACE(1),FSPACE(1),Y_SUPER,X,S_SUPER,S_STAR_SUPER,&
           Z,SK_DISCRETE_SUPER,SK_INTERP_SUPER)

       ELSE


        CALL APPROX (I, X, Z, A, FSPACE(IS6), FSPACE(1), ISPACE(1), FSPACE(IS5),&
                    FSPACE(IS4), ISPACE(2), ISPACE(3), ISPACE(5), ISPACE(8), & 
                    ISPACE(4), 2, DUMMY, 0)


       ENDIF

  
      RETURN


      END SUBROUTINE APPSLN


      SUBROUTINE APPROX (I, X, ZVAL, A, COEF, XI, N, Z, DMZ, K, &
                         NCOMP, MMAX, M, MSTAR, MODE, DMVAL, MODM )
!
!**********************************************************************
!
!   purpose
!                                    (1)       (m1-1)     (mncomp-1)
!           evaluate z(u(x))=(u (x),u (x),...,u  (x),...,u  (x)      )
!                              1     1         1          mncomp
!           at one point x.
!
!   variables
!     a      - array of mesh independent rk-basis coefficients
!     basm   - array of mesh dependent monomial coefficients
!     xi     - the current mesh (having n subintervals)
!     z      - the current solution vector
!     dmz    - the array of mj-th derivatives of the current solution
!     mode   - determines the amount of initialization needed
!            = 4  forms z(u(x)) using z, dmz and ha
!            = 3  as in =4, but computes local rk-basis
!            = 2  as in =3, but determines i such that
!                       xi(i) .le. x .lt. xi(i+1) (unless x=xi(n+1))
!            = 1  retrieve  z=z(u(x(i)))  directly
!
!**********************************************************************
!
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION ZVAL(1), DMVAL(1), XI(1), M(1), A(7,1), DM(7)
      DIMENSION Z(1), DMZ(1), BM(4), COEF(1)
!
      COMMON /COLOUT/ PRECIS, IOUT, IPRINT
!
!------------------------------------------------------------------------------------
!     Common Block
!
!     colout- 1) Precis is a constant = 100*macine unit roundoff
!             2) Iout is a constant = 6 for writing output to the screen
!             3) Iprint is equal to ipar(7) and therefore iprint
!         
!                = -1 for full diagnostic printout
!                = 0 for selected printout
!                = 1 for no printout
!      
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Appsln in Main Program or Colnew Wrapper
!               
!     Calls to : Rkbas
!
!-----------------------------------------------------------------------------------------------------------------------
      GO TO (10, 30, 80, 90), MODE
!
!...  mode = 1 , retrieve  z( u(x) )  directly for x = xi(i).
!
   10 X  = XI(I)
      IZ = (I-1) * MSTAR
      DO 20 J = 1, MSTAR
        IZ = IZ + 1
        ZVAL(J) = Z(IZ)
   20 CONTINUE
      RETURN
!
!...  mode = 2 ,  locate i so  xi(i) .le. x .lt. xi(i+1)
!
   30 CONTINUE
      IF ( X .GE. XI(1)-PRECIS .AND. X .LE. XI(N+1)+PRECIS ) &
                                                    GO TO 40
      IF (IPRINT .LT. 1)  WRITE(IOUT,900) X, XI(1), XI(N+1)
      IF ( X .LT. XI(1) )  X = XI(1)
      IF ( X .GT. XI(N+1) )  X = XI(N+1)
   40 IF ( I .GT. N .OR. I .LT. 1 )  I = (N+1) / 2
      ILEFT = I
      IF ( X .LT. XI(ILEFT) )                       GO TO 60
      DO 50 L = ILEFT, N
           I = L
           IF ( X .LT. XI(L+1) )                    GO TO 80
   50 CONTINUE
      GO TO 80
   60 IRIGHT = ILEFT - 1
      DO 70 L = 1, IRIGHT
           I = IRIGHT + 1 - L
           IF ( X .GE. XI(I) )                      GO TO 80
   70 CONTINUE
!
!...  mode = 2 or 3 , compute mesh independent rk-basis.
!
   80 CONTINUE
      S = (X - XI(I)) / (XI(I+1) - XI(I))
      CALL RKBAS ( S, COEF, K, MMAX, A, DM, MODM )
!
!...  mode = 2, 3, or 4 , compute mesh dependent rk-basis.
!
   90 CONTINUE
      BM(1) = X - XI(I)
      DO 95 L = 2, MMAX
         BM(L) = BM(1) / DFLOAT(L)
   95 CONTINUE
!
!...  evaluate  z( u(x) ).
!
  100 IR = 1
      IZ = (I-1) * MSTAR + 1
      IDMZ = (I-1) * K * NCOMP
      DO 140 JCOMP = 1, NCOMP
          MJ = M(JCOMP)
          IR = IR + MJ
          IZ = IZ + MJ
          DO 130 L = 1, MJ
             IND = IDMZ + JCOMP
             ZSUM = 0.D0
             DO 110 J = 1, K
               ZSUM = ZSUM  +  A(J,L) * DMZ(IND)
  110        IND = IND + NCOMP
             DO 120 LL = 1, L
               LB = L + 1 - LL
  120          ZSUM = ZSUM * BM(LB)  +  Z(IZ-LL)
  130     ZVAL(IR-L) = ZSUM
  140 CONTINUE
      IF ( MODM .EQ. 0 )                            RETURN
!
!...  for modm = 1 evaluate  dmval(j) = mj-th derivative of uj.
!
      DO 150 JCOMP = 1, NCOMP
  150 DMVAL(JCOMP) = 0.D0
      IDMZ = IDMZ + 1
      DO 170 J = 1, K
         FACT = DM(J)
         DO 160 JCOMP = 1, NCOMP
            DMVAL(JCOMP) = DMVAL(JCOMP)  +  FACT * DMZ(IDMZ)
            IDMZ = IDMZ + 1
  160    CONTINUE
  170 CONTINUE
      RETURN
!--------------------------------------------------------------------
  900 FORMAT(37H ****** DOMAIN ERROR IN APPROX ****** &
             /4H X =,D20.10, 10H   ALEFT =,D20.10, &
             11H   ARIGHT =,D20.10)


      END SUBROUTINE APPROX




      SUBROUTINE RKBAS (S, COEF, K, M, RKB, DM, MODE)
!
!**********************************************************************
!
!   purpose
!           evaluate mesh independent runge-kutta basis for given s
!
!   variables
!     s      - argument, i.e. the relative position for which
!              the basis is to be evaluated ( 0. .le. s .le. 1. ).
!     coef   - precomputed derivatives of the basis
!     k      - number of collocatin points per subinterval
!     m      - maximal order of the differential equation
!     rkb    - the runge-kutta basis (0-th to (m-1)-th derivatives )
!     dm     - basis elements for m-th derivative
!
!**********************************************************************
!
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Consts, Approx
!               
!     Calls to : No subroutines
!
!-----------------------------------------------------------------------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION COEF(K,1), RKB(7,1), DM(1), T(10)
!
      IF ( K .EQ. 1 )                            GO TO 70
      KPM1 = K + M - 1
      DO 10 I = 1, KPM1
   10   T(I) = S / DFLOAT(I)
      DO 40 L = 1, M
         LB = K + L + 1
         DO 30 I = 1, K
           P = COEF(1,I)
           DO 20 J = 2, K
             P = P * T(LB-J)  + COEF(J,I)
   20      CONTINUE
           RKB(I,L) = P
   30    CONTINUE
   40 CONTINUE
      IF ( MODE .EQ. 0 )                         RETURN
      DO 60 I = 1, K
         P = COEF(1,I)
         DO 50 J = 2, K
   50       P = P * T(K+1-J) + COEF(J,I)
         DM(I) = P
   60 CONTINUE
      RETURN
   70 RKB(1,1) = 1.0D0
      DM(1) = 1.0D0
      RETURN
    

      END SUBROUTINE RKBAS



      SUBROUTINE VMONDE ( RHO, COEF, K )
!
!**********************************************************************
!
!   purpose
!          solve vandermonde system v * x = e
!          with  v(i,j) = rho(j)**(i-1)/(i-1)! .
!
!**********************************************************************
!
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Consts
!               
!     Calls to : No subroutines
!
!-----------------------------------------------------------------------------------------------------------------------
      INTEGER K, I,IFAC,J,KM1,KMI
      DOUBLE PRECISION RHO(K), COEF(K)
!
      IF ( K .EQ. 1 )                             RETURN
      KM1 = K - 1
      DO 10 I = 1, KM1
         KMI = K - I
         DO 10 J = 1, KMI
           COEF(J) = (COEF(J+1) - COEF(J)) / (RHO(J+I) - RHO(J))
  10  CONTINUE
!
      IFAC = 1
      DO 40 I = 1, KM1
         KMI = K + 1 - I
         DO 30 J = 2, KMI
  30        COEF(J) = COEF(J) - RHO(J+I-1) * COEF(J-1)
         COEF(KMI) = DFLOAT(IFAC) * COEF(KMI)
         IFAC = IFAC * I
  40  CONTINUE
      COEF(1) = DFLOAT(IFAC) * COEF(1)
      RETURN
     

      END SUBROUTINE VMONDE




      SUBROUTINE HORDER (I, UHIGH, HI, DMZ, NCOMP, K)
!
!**********************************************************************
!
!   purpose
!           determine highest order (piecewise constant) derivatives
!           of the current collocation solution
!
!   variables
!     hi     - the stepsize, hi = xi(i+1) - xi(i)
!     dmz    - vector of mj-th derivative of the solution
!     uhigh  - the array of highest order (piecewise constant)
!              derivatives of the approximate solution on
!              (xi(i),xi(i+1)), viz,
!                          (k+mj-1)
!              uhigh(j) = u   (x)    on (xi(i),xi(i+1))
!                          j
!
!**********************************************************************
!
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION UHIGH(1), DMZ(1)
!
      COMMON /COLLOC/ RHO(7), COEF(49)
!
!-----------------------------------------------------------------------------------------------------------------
!     Common Block
!
!     colloc - 1) Rho are the k collocation points on (0,1)
!              2) Coef are the precomputed derivatives of the basis
!
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Newmsh
!               
!     Calls to : No subroutines
!
!-----------------------------------------------------------------------------------------------------------------------

      DN = 1.D0 / HI**(K-1)
!
!...  loop over the ncomp solution components
!
      DO 10 ID = 1, NCOMP
         UHIGH(ID) = 0.D0
   10 CONTINUE
      KIN = 1
      IDMZ = (I-1) * K * NCOMP + 1
      DO 30 J = 1, K
         FACT = DN * COEF(KIN)
         DO 20 ID = 1, NCOMP
            UHIGH(ID) = UHIGH(ID)  +  FACT * DMZ(IDMZ)
            IDMZ = IDMZ + 1
   20    CONTINUE
         KIN = KIN + K
   30 CONTINUE
      RETURN


      END SUBROUTINE HORDER



      SUBROUTINE DMZSOL (KD, MSTAR, N, V, Z, DMZ)
!
!**********************************************************************
!
!   purpose
!          compute dmz in a blockwise manner
!          dmz(i) = dmz(i)  +  v(i) * z(i), i = 1,...,n
!
!**********************************************************************
!
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Lsyslv
!               
!     Calls to : No subroutines
!
!-----------------------------------------------------------------------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION V(KD,1), DMZ(KD,1), Z(1)
!
      JZ = 1
      DO 30 I = 1, N
         DO 20 J = 1, MSTAR
            FACT = Z(JZ)
            DO 10 L = 1, KD
               DMZ(L,I) = DMZ(L,I)  +  FACT * V(L,JZ)
   10       CONTINUE
            JZ = JZ + 1
   20    CONTINUE
   30 CONTINUE
      RETURN


      END SUBROUTINE DMZSOL



!----------------------------------------------------------------------
!                            p a r t  5
!          we list here a modified (column oriented, faster)
!          version of the package solveblok of de boor - weiss [5].
!          we also give a listing of the linpack
!          routines dgefa und dgesl used by colnew.
!----------------------------------------------------------------------
!
      SUBROUTINE FCBLOK (BLOKS, INTEGS, NBLOKS, IPIVOT, SCRTCH, INFO)
!
!
!     calls subroutines  factrb  and  shiftb .
!
!     fcblok  supervises the plu factorization with pivoting of
!     scaled rows of the almost block diagonal matrix stored in the
!     arrays  bloks  and  integs .
!
!     factrb = subprogram which carries out steps 1,...,last of gauss
!            elimination (with pivoting) for an individual block.
!     shiftb = subprogram which shifts the remaining rows to the top of
!            the next block
!
!     parameters
!      bloks   an array that initially contains the almost block diago-
!            nal matrix  a  to be factored, and on return contains the
!            computed factorization of  a .
!      integs  an integer array describing the block structure of  a .
!      nbloks  the number of blocks in  a .
!      ipivot  an integer array of dimension   sum (integs(3,n) ; n=1,
!            ...,nbloks) which, on return, contains the pivoting stra-
!            tegy used.
!      scrtch  work area required, of length  max (integs(1,n) ; n=1,
!            ...,nbloks).
!      info    output parameter;
!            = 0  in case matrix was found to be nonsingular.
!            otherwise,
!            = n if the pivot element in the nth gauss step is zero.
!
!**********************************************************************
!
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Lsyslv
!               
!     Calls to : Shiftb, Factrb
!
!-----------------------------------------------------------------------------------------------------------------------
      INTEGER INTEGS(3,NBLOKS),IPIVOT(1),INFO, I,INDEX,INDEXN,LAST,&
              NCOL,NROW
      DOUBLE PRECISION BLOKS(1),SCRTCH(1)
      INFO = 0
      INDEXX = 1
      INDEXN = 1
      I = 1
!
!...  loop over the blocks.  i  is loop index
!
   10      INDEX = INDEXN
           NROW = INTEGS(1,I)
           NCOL = INTEGS(2,I)
           LAST = INTEGS(3,I)
!
!...       carry out elimination on the i-th block until next block
!...       enters, i.e., for columns 1,...,last  of i-th block.
!
           CALL FACTRB ( BLOKS(INDEX), IPIVOT(INDEXX), SCRTCH, NROW,&
                         NCOL, LAST, INFO)
!
!...       check for having reached a singular block or the last block
!
           IF ( INFO .NE. 0 )                       GO TO 20
           IF ( I .EQ. NBLOKS )                     RETURN
           I = I+1
           INDEXN = NROW * NCOL + INDEX
           INDEXX = INDEXX + LAST
!
!...       put the rest of the i-th block onto the next block
!
           CALL SHIFTB ( BLOKS(INDEX), NROW, NCOL, LAST, &
                         BLOKS(INDEXN), INTEGS(1,I), INTEGS(2,I) )
      GO TO 10
   20 INFO = INFO + INDEXX - 1
      RETURN


      END SUBROUTINE FCBLOK



      SUBROUTINE FACTRB ( W, IPIVOT, D, NROW, NCOL, LAST, INFO)
!
!********************************************************************
!
!     adapted from p.132 of  element.numer.analysis  by conte-de boor
!
!     constructs a partial plu factorization, corresponding to steps
!      1,..., last   in gauss elimination, for the matrix  w  of
!      order ( nrow ,  ncol ), using pivoting of scaled rows.
!
!     parameters
!       w       contains the (nrow,ncol) matrix to be partially factored
!               on input, and the partial factorization on output.
!       ipivot  an integer array of length last containing a record of
!               the pivoting strategy used; explicit interchanges
!               are used for pivoting.
!       d       a work array of length nrow used to store row sizes
!               temporarily.
!       nrow    number of rows of w.
!       ncol    number of columns of w.
!       last    number of elimination steps to be carried out.
!       info    on output, zero if the matrix is found to be non-
!               singular, in case a zero pivot was encountered in row
!               n,  info = n on output.
!
!**********************************************************************
!
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Fcblok
!               
!     Calls to : No subroutines
!
!-----------------------------------------------------------------------------------------------------------------------
      INTEGER IPIVOT(NROW),NCOL,LAST,INFO, I,J,K,L,KP1
      DOUBLE PRECISION W(NROW,NCOL),D(NROW), COLMAX,T,S
      DOUBLE PRECISION DABS,DMAX1
!
!...  initialize  d
!
      DO 10 I = 1, NROW
        D(I) = 0.D0
   10 CONTINUE
      DO 20 J = 1, NCOL
        DO 20 I = 1, NROW
          D(I) = DMAX1( D(I) , DABS(W(I,J)))
   20 CONTINUE
!
!...  gauss elimination with pivoting of scaled rows, loop over
!...  k=1,.,last
!
      K = 1
!
!...  as pivot row for k-th step, pick among the rows not yet used,
!...  i.e., from rows  k ,..., nrow , the one whose k-th entry
!...  (compared to the row size) is largest. then, if this row
!...  does not turn out to be row k, interchange row k with this
!...  particular row and redefine ipivot(k).
!
   30      CONTINUE
           IF ( D(K) .EQ. 0.D0 )                    GO TO 90
           IF (K .EQ. NROW)                         GO TO 80
           L = K
           KP1 = K+1
           COLMAX = DABS(W(K,K)) / D(K)
!
!...       find the (relatively) largest pivot
!
           DO 40 I = KP1, NROW
             IF ( DABS(W(I,K)) .LE. COLMAX * D(I) ) GO TO 40
             COLMAX = DABS(W(I,K)) / D(I)
             L = I
   40      CONTINUE
           IPIVOT(K) = L
           T = W(L,K)
           S = D(L)
           IF ( L .EQ. K )                          GO TO 50
             W(L,K) = W(K,K)
             W(K,K) = T
             D(L) = D(K)
             D(K) = S
   50      CONTINUE
!
!...       if pivot element is too small in absolute value, declare
!...       matrix to be noninvertible and quit.
!
           IF ( DABS(T)+D(K) .LE. D(K) )            GO TO 90
!
!...       otherwise, subtract the appropriate multiple of the pivot
!...       row from remaining rows, i.e., the rows (k+1),..., (nrow)
!...       to make k-th entry zero. save the multiplier in its place.
!...       for high performance do this operations column oriented.
!
           T = -1.0D0/T
           DO 60 I = KP1, NROW
   60        W(I,K) = W(I,K) * T
           DO 70 J=KP1,NCOL
             T = W(L,J)
             IF ( L .EQ. K )                        GO TO 62
               W(L,J) = W(K,J)
               W(K,J) = T
   62        IF ( T .EQ. 0.D0 )                     GO TO 70
             DO 64 I = KP1, NROW
   64           W(I,J) = W(I,J) + W(I,K) * T
   70      CONTINUE
           K = KP1
!
!...       check for having reached the next block.
!
           IF ( K .LE. LAST )                       GO TO 30
      RETURN
!
!...  if  last  .eq. nrow , check now that pivot element in last row
!...  is nonzero.
!
   80 IF( DABS(W(NROW,NROW))+D(NROW) .GT. D(NROW) ) RETURN
!
!...  singularity flag set
!
   90 INFO = K
      RETURN


      END SUBROUTINE FACTRB



      SUBROUTINE SHIFTB (AI, NROWI, NCOLI, LAST, AI1, NROWI1, NCOLI1)
!
!*********************************************************************
!
!     shifts the rows in current block, ai, not used as pivot rows, if
!     any, i.e., rows  (last+1),..., (nrowi), onto the first mmax =
!      = nrow-last  rows of the next block, ai1, with column last+j of
!      ai  going to column j , j=1,...,jmax=ncoli-last. the remaining
!     columns of these rows of ai1 are zeroed out.
!
!                                picture
!
!          original situation after         results in a new block i+1
!          last = 2 columns have been       created and ready to be
!          done in factrb (assuming no      factored by next factrb
!          interchanges of rows)            call.
!                      1
!                 x  x 1x  x  x           x  x  x  x  x
!                      1
!                 0  x 1x  x  x           0  x  x  x  x
!     block i          1                       ---------------
!     nrowi = 4   0  0 1x  x  x           0  0 1x  x  x  0  01
!     ncoli = 5        1                       1             1
!     last = 2    0  0 1x  x  x           0  0 1x  x  x  0  01
!     -------------------------------          1             1   new
!                      1x  x  x  x  x          1x  x  x  x  x1  block
!                      1                       1             1   i+1
!     block i+1        1x  x  x  x  x          1x  x  x  x  x1
!     nrowi1= 5        1                       1             1
!     ncoli1= 5        1x  x  x  x  x          1x  x  x  x  x1
!     -------------------------------          1-------------1
!                      1
!
!*********************************************************************
!
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Fcblok
!               
!     Calls to : No subroutines
!
!-----------------------------------------------------------------------------------------------------------------------
      INTEGER LAST, J,JMAX,JMAXP1,M,MMAX
      DOUBLE PRECISION AI(NROWI,NCOLI),AI1(NROWI1,NCOLI1)
      MMAX = NROWI - LAST
      JMAX = NCOLI - LAST
      IF (MMAX .LT. 1 .OR. JMAX .LT. 1)             RETURN
!
!...  put the remainder of block i into ai1
!
      DO 10 J=1,JMAX
           DO 10 M=1,MMAX
   10 AI1(M,J) = AI(LAST+M,LAST+J)
      IF (JMAX .EQ. NCOLI1)                         RETURN
!
!...  zero out the upper right corner of ai1
!
      JMAXP1 = JMAX + 1
      DO 20 J=JMAXP1,NCOLI1
           DO 20 M=1,MMAX
   20 AI1(M,J) = 0.D0
      RETURN

      END SUBROUTINE SHIFTB



      SUBROUTINE SBBLOK ( BLOKS, INTEGS, NBLOKS, IPIVOT, X )
!
!**********************************************************************
!
!     calls subroutines  subfor  and  subbak .
!
!     supervises the solution (by forward and backward substitution) of
!     the linear system  a*x = b  for x, with the plu factorization of
!     a  already generated in  fcblok .  individual blocks of
!     equations are solved via  subfor  and  subbak .
!
!    parameters
!       bloks, integs, nbloks, ipivot    are as on return from fcblok.
!       x       on input: the right hand side, in dense storage
!               on output: the solution vector
!
!*********************************************************************
!
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Lsyslv
!               
!     Calls to : Subfor, Subbak
!
!-----------------------------------------------------------------------------------------------------------------------
      INTEGER INTEGS(3,NBLOKS),IPIVOT(1), I,INDEX,INDEXX,J,LAST, &
              NBP1,NCOL,NROW
      DOUBLE PRECISION BLOKS(1), X(1)
!
!...  forward substitution pass
!
      INDEX = 1
      INDEXX = 1
      DO 10 I = 1, NBLOKS
           NROW = INTEGS(1,I)
           LAST = INTEGS(3,I)
           CALL SUBFOR ( BLOKS(INDEX), IPIVOT(INDEXX), NROW, LAST,&
                         X(INDEXX) )
           INDEX = NROW * INTEGS(2,I) + INDEX
   10 INDEXX = INDEXX + LAST
!
!...  back substitution pass
!
      NBP1 = NBLOKS + 1
      DO 20 J = 1, NBLOKS
           I = NBP1 - J
           NROW = INTEGS(1,I)
           NCOL = INTEGS(2,I)
           LAST = INTEGS(3,I)
           INDEX = INDEX - NROW * NCOL
           INDEXX = INDEXX - LAST
   20 CALL SUBBAK ( BLOKS(INDEX), NROW, NCOL, LAST, X(INDEXX) )
      RETURN


      END SUBROUTINE SBBLOK



      SUBROUTINE SUBFOR ( W, IPIVOT, NROW, LAST, X )
!
!**********************************************************************
!
!     carries out the forward pass of substitution for the current
!     block, i.e., the action on the right side corresponding to the
!     elimination carried out in  factrb  for this block.
!
!    parameters
!       w, ipivot, nrow, last  are as on return from factrb.
!       x(j)  is expected to contain, on input, the right side of j-th
!             equation for this block, j=1,...,nrow.
!       x(j)  contains, on output, the appropriately modified right
!             side of equation (j) in this block, j=1,...,last and
!             for j=last+1,...,nrow.
!
!*********************************************************************
!
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Sbblok
!               
!     Calls to : No subroutines
!
!-----------------------------------------------------------------------------------------------------------------------
      INTEGER IPIVOT(LAST), IP,K,KP1,LSTEP
      DOUBLE PRECISION W(NROW,LAST), X(NROW), T
!
      IF ( NROW .EQ. 1 )                            RETURN
      LSTEP = MIN0( NROW-1 , LAST )
      DO 20 K = 1, LSTEP
           KP1 = K + 1
           IP = IPIVOT(K)
           T = X(IP)
           X(IP) = X(K)
           X(K) = T
           IF ( T .EQ. 0.D0 )                       GO TO 20
           DO 10 I = KP1, NROW
   10         X(I) = X(I) + W(I,K) * T
   20 CONTINUE
   30 RETURN


      END SUBROUTINE SUBFOR


      SUBROUTINE SUBBAK ( W, NROW, NCOL, LAST, X )
!
!*********************************************************************
!
!     carries out backsubstitution for current block.
!
!    parameters
!       w, ipivot, nrow, ncol, last  are as on return from factrb.
!       x(1),...,x(ncol)  contains, on input, the right side for the
!               equations in this block after backsubstitution has been
!               carried up to but not including equation (last).
!               means that x(j) contains the right side of equation (j)
!               as modified during elimination, j=1,...,last, while
!               for j .gt. last, x(j) is already a component of the
!               solution vector.
!       x(1),...,x(ncol) contains, on output, the components of the
!               solution corresponding to the present block.
!
!*********************************************************************
!
!----------------------------------------------------------------------------------------------------------------------
!
!     Calls by : Sbblok
!               
!     Calls to : No subroutines
!
!-----------------------------------------------------------------------------------------------------------------------
      INTEGER  LAST,  I,J,K,KM1,LM1,LP1
      DOUBLE PRECISION W(NROW,NCOL),X(NCOL), T
!
      LP1 = LAST + 1
      IF ( LP1 .GT. NCOL )                          GO TO 30
      DO 20 J = LP1, NCOL
         T = - X(J)
         IF ( T .EQ. 0.D0 )                         GO TO 20
         DO 10 I = 1, LAST
   10       X(I) = X(I) + W(I,J) * T
   20 CONTINUE
   30 IF ( LAST .EQ. 1 )                            GO TO 60
      LM1 = LAST - 1
      DO 50 KB = 1, LM1
        KM1 = LAST - KB
        K = KM1 + 1
        X(K) = X(K)/W(K,K)
        T = - X(K)
        IF ( T .EQ. 0.D0 )                          GO TO 50
        DO 40 I = 1, KM1
   40     X(I) = X(I) + W(I,K) * T
   50 CONTINUE
   60 X(1) = X(1)/W(1,1)
      RETURN

      END SUBROUTINE SUBBAK


      SUBROUTINE DFDY_DIVDIFF(X,Z,NODE,NSOL,FSUB,DF)

        ! Compute the divided difference approximation to d FSUB/dy. FSUB is the user routine which
        ! evaluates the right hand side of the differential equation.
        ! The derivatives are to be evaluated at a mesh point X and the corresponding solution
        ! approximation Z. NODE is the number of differential equations. NSOL are the number of solution
        ! components in Z. The derivative values are returned in DF
        !---------------------------------------------------------------------------------------------------
        ! Epsilon function is used to determine the round off error of a double precision number       
        INTRINSIC EPSILON

        ! Input arguments
        DOUBLE PRECISION :: X,Z(NSOL),DF(NODE,NSOL)

        ! Local arguments
        INTEGER::I,J !Loop index
        DOUBLE PRECISION,PARAMETER :: UROUND=EPSILON(1.0D0) !Double precision round off error
        DOUBLE PRECISION :: SAVEZI ! Ith component of Z
        DOUBLE PRECISION :: DELZI ! Denominator for divided difference
        DOUBLE PRECISION :: SQRTU ! Square root of the round off error
        DOUBLE PRECISION :: FZ(NODE) ! Value of FSUB at Z
        DOUBLE PRECISION :: FDELZ(NODE) ! Value of FSUB at Z+DELZI

        EXTERNAL FSUB

        !----------------------------------------------------------------------------------------------------
        SQRTU=SQRT(UROUND)

        CALL FSUB(X,Z,FZ)

        DO I=1,NSOL ! Start of loop I=1,NSOL

           SAVEZI=Z(I)

          ! Determine the appropriate perturbation of Z(I) and perturb Z(I) by this amount.
          IF(ABS(SAVEZI)> 0.0D0)THEN
            DELZI=SQRTU*ABS(SAVEZI)
          ELSE
            DELZI=SQRTU
          ENDIF

          Z(I)=SAVEZI+DELZI

          ! Evaluate FSUB for input X and perturbed Z
          CALL FSUB(X,Z,FDELZ)

          ! Compute the divided dfference. Store in DF.
          DO J=1,NODE

            DF(J,I)=(FDELZ(J)-FZ(J))/DELZI

          ENDDO

          ! Restore Z(I) value
          Z(I)=SAVEZI

        ENDDO ! End of loop I=1,NSOL

       END SUBROUTINE DFDY_DIVDIFF 

       SUBROUTINE DBCDY_DIVDIFF(I,Z,NSOL,GSUB,DG)

        ! Compute divided difference approximation to d GSUB/DY. GSUB is the user routine which
        ! evaluates the boundary conditions. The derivative is evaluated at the ITH boundary condition
        ! given the solution approximation Z at the endpoint. NSOL is the number of solution components in
        ! Z.The partial derivatives are returned in DG.
        !-------------------------------------------------------------------------------------------------
        INTRINSIC EPSILON

        ! Input Arguments
        INTEGER::I
        DOUBLE PRECISION :: Z(NSOL)

        ! Output Arguments
        DOUBLE PRECISION :: DG(NSOL)
        
        ! Local Variables
        INTEGER::K,J  ! Loop index
        DOUBLE PRECISION,PARAMETER :: UROUND=EPSILON(1.0D0)
        DOUBLE PRECISION :: SAVEZK ! Kth component of Z. 
        DOUBLE PRECISION :: DELZK ! Denominator fro divided difference
        DOUBLE PRECISION :: SQRTU ! Square root of round off error
        DOUBLE PRECISION :: GZ ! Value of BCSUB at Z
        DOUBLE PRECISION :: GDELZ ! Value of BCSUB at Z + DELZK

        EXTERNAL GSUB

        !------------------------------------------------------------------------------------------------

        SQRTU=SQRT(UROUND)

        CALL GSUB(I,Z,GZ)

        DO K=1,NSOL

         SAVEZK=Z(K)

         ! Determine appropriate perturbation of Z(K)
         IF(ABS(SAVEZK)> 0.0D0)THEN
            DELZK=SQRTU*ABS(SAVEZK)
         ELSE
            DELZK=SQRTU
         ENDIF

         Z(K)=SAVEZK+DELZK

         ! Evaluate GSUB for input I and perturbed Z. GDELZ is returned.
         CALL GSUB(I,Z,GDELZ)

         DG(K)=(GDELZ-GZ)/DELZK

         ! Restore Z(K)
         Z(K)=SAVEZK

        ENDDO

       END SUBROUTINE DBCDY_DIVDIFF 



! ***Beginning of SCSA technology*************************************************




       SUBROUTINE EVALUATE(NEQNS,NSUB,MESH,Y,T,S,S_STAR,Z, &
                               K_DISCRETE,K_INTERP)
!
!  ***OVERVIEW***
!                
!          This routine uses the current solution approximation, 'Y',
!       obtained on the mesh stored in 'mesh', plus the stages of the 
!       underlying discrete RK method plus some new stages associated with
!       the interpolant in order to evaluate the interpolant.  This interpolant 
!       is evaluated at 't', with the value returned in 'z', and the value of 
!       the first derivative returned in 'z_prime'.
!
!          The interpolant is computed as follows: Given 't', assumed to 
!       be somewhere in the problem interval, the subinterval of 'mesh' that 
!       contains 't' is determined first, i.e. find i such that mesh(i-1) 
!       =< t < mesh(i) for t < mesh(Nsub) or i = Nsub for t = mesh(Nsub).
!       The fraction of the way 't' is through the ith subinterval, a quantity
!       which we call 'tau',  is also determined during this step. Then the 
!       interpolant on this ith subinterval is evaluated - the interpolant has
!       the form z(mesh(i-1) + tau*hi) = yim1 + hi sum br(tau)*kr. The weights,
!       br(tau), can be computed with a call to 'interp_weights'. 'Nsub' is 
!       the number of subintervals.
!
!          The remaining information required to define the interpolant 
!       on the ith subinterval is the current solution approximation at the 
!       (i-1)st mesh point, 'yim1', which can be obtained from 'Y', and the
!       stages associated with the ith subinterval. The first 's' stages,
!       contained in 'k_discrete', come from the discrete formula and
!       remaining 's_star - s' stages, contained in 'k_interp', are necessary
!       for the interpolant. Each stage is a vector of length 'neqns' where 
!       'neqns' is the number of differential equations. 
!
!          'k_discrete' is the vector containing the discrete stages for
!       all subintervals. The ith set of 's*neqns' locations of 'k_discrete' 
!       contains the 's' stages, each of length 'neqns'corresponding to the 
!       ith subinterval. 'k_interp' plays a similar role for the 's_star-s'
!       stages associated with the interpolant.
!

!------------------------------------------------------------------------------
!
!  ***declaration of constants***
!
        INTEGER MXS 
        PARAMETER(MXS=10)
!
!       'Mxs'    is the maximum value for the total number of stages
!                of the interpolant.
!
!  ***declaration of parameters***
!
!      imports:
!
        INTEGER NEQNS
        INTEGER NSUB
        DOUBLE PRECISION MESH(0:NSUB)
        DOUBLE PRECISION Y(NEQNS*(NSUB+1))
        DOUBLE PRECISION T
!
!       'neqns' is the number of differential equations.
!
!       'Nsub'  is the number of subintervals
!
!       'mesh'  is the list of points defining the subintervals
!
!       'Y'     contains the discrete solution at the mesh points
!               the ith set of 'neqns' locations of 'Y' contains 
!               the solution approximation at the ith mesh point.
!
!       't'     the point at which the interpolant is to be evaluated.
!
        DOUBLE PRECISION K_DISCRETE(MXS*NEQNS*NSUB)
!                
!       'k_discrete' contains the 's' discrete stage values for each
!                               subinterval.
!
        DOUBLE PRECISION K_INTERP(MXS*NEQNS*NSUB)
!
!       'k_interp' contains the 's_star-s' extra stage values for each
!                  subinterval need to form the interpolant.
!
!      exports:
!
        DOUBLE PRECISION Z(NEQNS)
     
!
!       'z'       is the value of the interpolant at x.
!
!       'z_prime' is the value of the first derivative of the
!                 interpolant at x.
!
!  ***declaration of local variables***
!
        INTEGER I
        DOUBLE PRECISION HI
        DOUBLE PRECISION TAU
        DOUBLE PRECISION WEIGHTS(MXS)
  
!
!       'i'             index for subinterval containing 't', i ranges
!                       from 1 to 'Nsub'
!
!       'hi'            the size of the ith subinterval.
!
!       'tau'           equal to (t-mesh(i-1))/hi
!
!       'weights'       contains the weight polynomials of the interpolant
!                       evaluated at 'tau'.
!
!       'weights_prime' contains the first derivatives of the weight 
!                       polynomials of the interpolant evaluated at 'tau'.
!
!  ***declaration of variables from /RK_s/ common block***
        INTEGER S
!       's'          is the number of discrete stages.


!  ***declaration of variables from /interp_s_star/ common block***
        INTEGER S_STAR
!       's_star'   is the 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.
!
!-----------------------------------------------------------------------------
!      called by:'col_mirk',    user's main program 
!    
!
!      calls to: 'interval', 'col_weights', 'sum_stages'
!------------------------------------------------------------------------------
!
!         Call 'interval' to determine the subinterval containing t.
!         The index of the subinterval is returned in 'i'.
!
!
          CALL INTERVAL(NSUB,MESH,T,I)
!                                                                        
!
!         Compute 'tau' for use in the evaluation of the weight polynomials.
!
          HI = MESH(I) - MESH(I-1)
          TAU = (T - MESH(I-1))/HI
!
!         Setup the weights, evaluated at 'tau'.
!
          CALL COL_WEIGHTS(S_STAR, TAU, WEIGHTS)
!
!         We need to evaluate the interpolant and its derivative at 'tau'.
!
!         This is done by taking weighted sums of the stages, using the weights
!         we have just computed, and the stages corresponding to the ith
!         subinterval, and calling the 'sum_stages' routine.
!

          CALL SUM_STAGES(NEQNS,HI,Y((I-1)*NEQNS+1),S, &
                          K_DISCRETE((I-1)*S*NEQNS+1),S_STAR, &
                          K_INTERP((I-1)*(S_STAR-S)*NEQNS+1), &
                          WEIGHTS, Z)

!


          RETURN
          END SUBROUTINE EVALUATE

      SUBROUTINE INTERVAL(NSUB, MESH, T, I)
!      
!  ***OVERVIEW***
!
!         This routine finds i such that 'mesh(i-1) =< t < mesh(i)' for
!         t < mesh(Nsub). The value of Nsub is returned when t = mesh(Nsub).
!         'Nsub' is the number of subintervals and 'mesh' is the list of points
!         defining the subintervals. 
!
!----------------------------------------------------------------------------
!
!  ***declaration of parameters***
!
!      imports:
!             
          INTEGER NSUB
          DOUBLE PRECISION MESH(0:NSUB), T                  
!
!         'Nsub'  is the number of subintervals.
!
!         'mesh'  is the list of points defining the subintervals.
!
!         't'     the point to be located within 'mesh'.
!
!         'meshhalfnew' the new mesh value resulting
!                       from halving the current interval.
!
!      exports:
!
          INTEGER I
!
!         'i'   t satisfies mesh(i-1) =< t < mesh(i)  for t < mesh(Nsub)
!               or i = Nsub for t = mesh(Nsub)
!
!  ***declaration of local variables***
!
          DOUBLE PRECISION     MESHHALFNEW
          INTEGER              INEW
          INTEGER              IOLD
          INTEGER              IRIGHT
          INTEGER              ILEFT
!
!         'meshhalfnew'  new mesh value resulting from halving
!                        the current interval
!
!         'inew'         current subinterval containing
!                        meshhalfnew
!
!         'iold'         previous subinterval containing 
!                        meshhalfnew
!
!         'ileft'        lefthand mesh point bounding
!                        meshhalfnew
!
!         'iright'       righthand mesh point bounding
!                        meshhalfnew
!
!--------------------------------------------------------------------------- 
!      called by: 'interp_eval' or 'evaluate'
!---------------------------------------------------------------------------
!
!     Handle special cases for the location of t first
!
      IF ((T .LT. MESH(0)) .OR. (T .GT. MESH(NSUB))) THEN
         WRITE(6,*)'ERROR-T IS OUTSIDE OF PROBLEM INTERVAL'
      ENDIF
!      
      IF (T .EQ. MESH(NSUB)) THEN
                  I = NSUB
                  RETURN
      END IF
!
      IF (T .EQ. MESH(0)) THEN
                  I = 1
                  RETURN
      END IF
!
!     At this point we know that mesh(0) < t < mesh(Nsub).
!
!     *********A binary search is used here***************
!             
      INEW = (NSUB/2)
      MESHHALFNEW = MESH(INEW)
!
!     ***determine if t is equal to meshhalfnew***
!
      IF (T .EQ. MESHHALFNEW) THEN
         I = INEW + 1
         RETURN
      ENDIF
!
!     ****t is not equal to the meshpoint******
!
      IF (T .LT. MESHHALFNEW)THEN
!
!     ****t is on the lefthand side of the interval****
!
           DO WHILE (T .LT. MESHHALFNEW)
              IOLD = INEW
              INEW = IOLD/2
              MESHHALFNEW = MESH(INEW)
           ENDDO         
!  
           IRIGHT = IOLD
           ILEFT = INEW
!
       ELSE
!
           DO WHILE (T .GT. MESHHALFNEW)
!
!      ****t is on the righthand side of the interval****
!
              IOLD = INEW
              INEW = NINT(REAL((IOLD + (NSUB - IOLD)/2)) + 0.5)
              MESHHALFNEW = MESH(INEW)                                   
           ENDDO
!
           ILEFT  = IOLD
           IRIGHT = INEW
!
       ENDIF 
!     
       I = IRIGHT
!
!      ***determine if t is equal to meshhalfnew***
!
       IF (T .EQ. MESHHALFNEW) THEN
         I = INEW + 1
         RETURN
       ENDIF
!       
!      ****t is not equal to the meshpoint******
!
!      ****t is bounded by ileft and iright****
!
       DO WHILE (ABS(IRIGHT - ILEFT) .NE. 1)
!           
           INEW = (ILEFT + IRIGHT)/2
!
           IF(MESH(INEW) .LT. T) THEN
                 ILEFT = INEW
           ELSEIF (MESH(INEW) .GT. T) THEN
                 IRIGHT = INEW
           ELSE
                 I = INEW + 1
                 RETURN
           ENDIF
!
       ENDDO
!
       I = IRIGHT
!
!      ****t is in the ith subinterval****
! 
       RETURN
       END SUBROUTINE INTERVAL
          
       SUBROUTINE SETUP(NODE,NEQNS,ORDERS,TIM1,HI,YIM1,YI,S,KI_DISCRETE,&
                        S_STAR,KI_INTERP,Y,FSUB)
!
!  ***OVERVIEW***
!
!         The purpose of this routine is to perform some preliminary 
!      calculations to set up stages for later use in the evaluation of the 
!      interpolant. The interpolant is defined on a subinterval by subinterval 
!      basis, as follows: Zi(mesh(i-1) + hi*tau) = y + hi sum br(tau)*kr, where
!      the stages, kr, depend on information from the ith subinterval. The y 
!      value is available from the discrete solution and the value of hi, the 
!      size of the subinterval, is easily computed. The weight polynomials,
!      br(tau), can be evaluated through a call to the interp_weights routine.
!      'yim1' and 'yi' are the solution approximations at the left and right
!      endpoints, respectively. The first 's' stages of each subinterval were
!      already obtained during the computation of the discrete solution, so 
!      only the last 's_star - s' stages need to be computed now.
!
!         This routine will construct the 's_star-s' extra stages needed by 
!      the interpolant for a single subinterval. This routine takes as input
!      the vector, 'ki_discrete', containing the 's' stages associated with
!      the discrete formula on the ith subinterval, and then uses the coupling
!      coefficients for the extra stages, available from the /interp_coeffs/
!      common block, to construct the 's_star' new stages. These are returned
!      in the vector 'ki_interp'. Each stage is a vector of length 'neqns', 
!      where 'neqns' is the number of differential equations. 'tim1' is the 
!      left endpoint of the ith subinterval.
!
!         'y' is a work array.
!-------------------------------------------------------------------------------
!
!  ***declaration of constants***
!
          INTEGER  MXS 
          PARAMETER (MXS=10)

          INTEGER  MAXNODE
          PARAMETER (MAXNODE=50)

!
!         'Mxs'    is the maximum value for the total number of stages
!                  of the interpolant.
!
!  ***declaration of parameters***
!
!     imports:
!
          INTEGER            NEQNS,NODE
          INTEGER            ORDERS(NODE)
          DOUBLE PRECISION   F(MAXNODE)
          DOUBLE PRECISION   TIM1
          DOUBLE PRECISION   HI       
          DOUBLE PRECISION   YIM1(NEQNS)
          DOUBLE PRECISION   YI(NEQNS)
          INTEGER            S
          DOUBLE PRECISION   KI_DISCRETE(S*NEQNS)
          INTEGER            S_STAR
!
!         'neqns' is the number of ordinary differential equations.
!
!         'tim1'  is left endpoint of the ith subinterval.
!
!         'hi'    is the size of the current subinterval.
!
!         'yim1'  is the solution approximation at the left endpoint
!                 of the current subinterval.
! 
!         'yi'    is the solution approximation at the right endpoint of
!                 the current subinterval.
!
!         's'     is the number of discrete stages.
!
!         'ki_discrete' is a vector of 's*neqns' components. The jth 
!                       component is of length 'neqns' and corresponds
!                       to the jth of the 's' stage values associated
!                       with the discrete formula on the ith subinterval.
!
!         's_star'      is the 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.
!
!     exports:
!
          DOUBLE PRECISION KI_INTERP((S_STAR-S)*NEQNS)
!
!         'ki_interp' has the same structure as 'ki_discrete'.
!                     The jth vector component corresponds to the jth 
!                     of the new stages needed to form the interpolant.
!
!     work space:

          DOUBLE PRECISION Y(NEQNS)
!
!         'y'    is used to hold the argument for the function evaluation
!                for each new stage as it is being computed.
!
!
!  ***user - supplied subroutines***
!
!         fsub - This routine defines f(t,y) for the first order system
!                of differential equations, y' = f(t,y).
!
!  ***declaration of local variables***
!
          INTEGER          R
          INTEGER          J

!         'r,j'  loop indexes from 1 to 's_star-s'.
!
!  ***declaration of variables from common block /interp_coeff/***
!
!     imports:
!      
          DOUBLE PRECISION C_STAR(MXS)
          DOUBLE PRECISION V_STAR(MXS)
          DOUBLE PRECISION X_STAR(MXS*MXS)
!
!         The first component of 'c_star' contains the abscissa
!         for the first new stage. The second component of 'c_star' contains
!         the abscissa for the second new stage, and so on.
!
!         The first component of 'v_star' contains the blending coefficient
!         for the first new stage. The second component of 'v_star' contains
!         the blending coefficient for the second new stage, and so on.
!
!         The Runge-Kutta coefficient matrix, X_star, is related to the data
!         structure, 'x_star' as follows (the matrix is stored columnwise in
!         the vector):
!
!               X_star(i,j) = x_star( (j-1)*s + i )
!          
!         'x_star' contains the coupling coefficients for each of the new
!         stages

!
          COMMON /INTERP_COEFF/ C_STAR,V_STAR,X_STAR
!
          EXTERNAL FSUB
!------------------------------------------------------------------------------
!      called by: 'col_mirk' 
!      calls to : 'fsub', 'dcopy', 'daxpy'
!------------------------------------------------------------------------------
!
!  ***CONSTRUCTION OF THE NEW STAGES***
!
      DO 20 R = 1 , S_STAR-S
!
!	   write(6,*)'Coefficients for extra stage',r

!          initialize 'y' to zero
!
           CALL DCOPY(NEQNS,0.0D0,0,Y,1)
!
!          Accumulate the argument for the rth new stage in 'y'.
!
!          The contributions from the discrete stages.
!
           DO 5 J = 1 , S
!
!                The contribution from the jth discrete stage.
                 CALL DAXPY(NEQNS,X_STAR((J-1)*(S_STAR-S)+R),&
                                KI_DISCRETE((J-1)*NEQNS+1),1,Y,1)
!                 write(6,*) x_star((j-1)*(s_star-s)+r)
 5         CONTINUE
!
!          The contributions from the previous new stages.
!
            DO 10 J = 1 , R-1
!
!                 The contribution from the jth new stage.
                  CALL DAXPY(NEQNS,X_STAR((J+S-1)*(S_STAR-S)+R),&
                                KI_INTERP((J-1)*NEQNS+1),1,Y,1)
!                  write(6,*) x_star((j+s-1)*(s_star-s)+r)
 10         CONTINUE
!
!           Multiply the sum of the stages by 'hi'.
            CALL DSCAL(NEQNS,HI,Y,1)
!
!           Add on the contribution from yim1.
            CALL DAXPY(NEQNS,(1-V_STAR(R)),YIM1,1,Y,1)
!
!           Add on the contribution from yi.
            CALL DAXPY(NEQNS,V_STAR(R),YI,1,Y,1)
!
!           All contributions to the argument for rth new stage have
!           been made. Simply evaluate 'fsub' at this argument to get the
!           new stage. The result is copied into the appropriate locations
!           of 'ki_interp'.
!

            CALL FSUB(TIM1+C_STAR(R)*HI,Y,   &
                      KI_INTERP((R-1)*NEQNS+1))

     



!           end of calculation for rth new stage
!
 20   CONTINUE
!
      RETURN
      END SUBROUTINE SETUP


      SUBROUTINE SUM_STAGES(NEQNS,HI,YIM1,S,KI_DISCRETE,&
                  S_STAR,KI_INTERP,WEIGHTS,Z)
!
!  ***OVERVIEW***
!
!          Given the stages and corresponding weights, this routine
!      computes the value of the interpolant, 'z', and its first derivative, 
!      'z_prime', as follows: 
!            z = yim1 + hi sum weights(r)*ki(r),
!            z_prime = sum weights_prime(r)*ki(r),
!      where 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 only. Each stage is a vector of length 'neqns' where
!      'neqns' is the number of differential equations. 'hi' is the size of
!      the ith subinterval and 'yim1'  is the discrete solution at the left
!      endpoint of the ith 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.
!
!------------------------------------------------------------------------------
!
!   ***declaration of parameters***
!
!      imports:
!
            INTEGER NEQNS
            DOUBLE PRECISION HI
            DOUBLE PRECISION YIM1(NEQNS)
            INTEGER S
            DOUBLE PRECISION KI_DISCRETE(S*NEQNS)
            INTEGER S_STAR
            DOUBLE PRECISION KI_INTERP((S_STAR-S)*NEQNS)
            DOUBLE PRECISION WEIGHTS(S_STAR)
          
!
!           'neqns'  is the number of differential equations.
!
!           'hi'     is the size of the current (ith) subinterval.
!
!           'yim1'   is the discrete solution at the left endpoint
!                    of the ith subinterval.
!
!           's'      is the number of discrete stages.
!
!           'ki_discrete'   the set of 's' discrete stages associated with 
!                           the ith subinterval. 
!
!           's_star'        is the total number of stages.
!
!           'ki_interp'     the set of 's_star-s' interpolation stages 
!                           associated with the ith subinterval. 
!
!           'weights'       the weights associated with the stages needed to
!                           evaluate the interpolant at the appropriate
!                           point. The weights have aleady been evaluated at
!                           that point.
!           
!           'weights_prime' the weights associated with the stages needed to
!                           form the first derivative of the interpolant
!                           at the appropriate point. The weights have already
!                           been evaluated at this point.
!
!      exports:
!
            DOUBLE PRECISION Z(NEQNS)
         
!
!           'z'       is the value of the interpolant.
!
!           'z_prime' is the value of the first derivative of the 
!                     interpolant.
!
!   ***declaration of local variables***
!
            INTEGER J
!
!           'j'       loop index from 1 to 's', or 1 to 's_star-s'
!
!---------------------------------------------------------------------------
!      Called by: defect_estimate, interp_eval
!      Calls to : dcopy, daxpy
!---------------------------------------------------------------------------
!
!          Initialize all vectors to be used in this computation to zero.
!
           CALL DCOPY(NEQNS,0.0D0,0,Z,1)
          
!
!          We can now sum the appropriately weighted stages to compute
!          the interpolant 'z' and the first derivative of the 
!          interpolant, 'z_prime'.
!
!          Accumulate weighted sum of discrete stages for 'z' and 'z_prime'.
!

           DO 10 J = 1, S
!
                CALL DAXPY(NEQNS,WEIGHTS(J), &
                           KI_DISCRETE((J-1)*NEQNS+1),1,Z,1)
         
 10        CONTINUE
!
!          Accumulate weighted sum of new interpolant stages for 'z' 
!          and 'z_prime'.
!
           DO 20 J = 1,(S_STAR-S)
!
            CALL DAXPY(NEQNS,WEIGHTS(S+J), &
                           KI_INTERP((J-1)*NEQNS+1),1,Z,1)


 20        CONTINUE

!
!          The calculation of 'z_prime' is now complete. 
!          Scale the weighted sum of stages by hi, and then add
!          on yim1. This will complete the calculation of 'z'. 
! 

           CALL DSCAL(NEQNS,HI,Z,1)
           CALL DAXPY(NEQNS,1.0D0,YIM1,1,Z,1)
!            

           RETURN
           END SUBROUTINE SUM_STAGES


      SUBROUTINE TABLEAU_STAR(METHOD)
! 
! ***OVERVIEW***
!
!        The main purpose of this routine is to define the parameter which
!      gives,'s_star', the total number of stages for the interpolant and the 
!      coefficients for the 's_star-s' extra stages needed by the interpolant.
!      Once defined by this routine, the coefficients are then
!      available from the /interp_coeff/ common block, while 's_star' is
!      available from the /interp_s_star/ common block. The 'order' of the
!      interpolant is available through the /interp_order/ common block.
!      The particular coefficients that get defined here depend on which
!      'method' is being defined. For 'method=2', no extra stages are required
!      so this routine is not called; for 'method=3', one extra stage is needed;
!      for 'method=4', three extra stages are needed.
! 
!            The final portion of formula dependent information for the
!      interpolant involves
!      the expressions for the weight polynomials, arising as part of the
!      definition of the interpolant for the discrete formula. This information
!      is provided through the 'col_weights' routine which evaluates all
!      the weight polynomials and their derivatives at a given point. This
!      routine can be called by other parts of the code whenever the weight 
!      polynomials and derivatives are required - i.e. whenever an evaluation 
!      of the interpolant is required.
!
!-----------------------------------------------------------------------------
!
!   ***declaration of constants***
!
            INTEGER MXS
            PARAMETER(MXS=10) 
!
!   ***declaration of parameters
!
            INTEGER METHOD

!            'method' identification of the collocation scheme that is
!            embedded with the CRK scheme.

!   ***declaration of variables from common block /interp_s_star/***

            INTEGER S_STAR

!           's_star'  is the 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.

            COMMON /INTERP_S_STAR/ S_STAR
!
!   ***declaration of variables from common block /interp_coeff/***
!
!         exports:
!      
            DOUBLE PRECISION C_STAR(MXS), V_STAR(MXS), X_STAR(MXS*MXS)
!                                                                        
!           The first component of 'c_star' contains the abscissa
!           for the first new stage. The second component of 'c_star' contains
!           the abscissa for the second new stage, and so on.
!
!           The first component of 'v_star' contains the blending coefficient
!           for the first new stage. The second component of 'v_star' contains
!           the blending coefficient for the second new stage, and so on.
!
!           The matrix, X_star, is related to  the data structure, 'x_star'
!           as follows (the matrix is stored columnwise in the vector):
!
!                 X_star(i,j) = x_star( (j-1)*s + i )
!          
!           'x_star' contains the coupling coefficients for each of the new
!           stages
!
            COMMON /INTERP_COEFF/ C_STAR,V_STAR,X_STAR
!
!------------------------------------------------------------------------------
!      called by: 'cmirk'
!      calls to: 'dcopy'
!------------------------------------------------------------------------------
!
!  As part of SCSA technology, the more efficient approach is to hard code the tableau constants.



       IF (METHOD .EQ. 2) S_STAR = 4

       IF (METHOD .EQ. 3) THEN
!
!        There is one new stage.
!
         S_STAR = 6
!
!         C_STAR(1) = 3.D0/5.D0 

         c_star(1)= 0.6d0
!
!         V_STAR(1) = 3.D0/5.D0 
!
         v_star(1)=0.6d0


!        Column 1 of x_star
!
!         X_STAR(1) = - 78.D0/ 3125.D0
!
!          x_star(1)=-0.2496000000000000d-1

!          x_star(1)=-2.496000000000000000000000000000000000000d-2



           x_star(1)=-0.2496D-1

!        
!        Column 2 of x_star
!
!         X_STAR(2) = 72.D0/3125.D0
!
!          x_star(2)=0.2304000000000000d-1

!          x_star(2)=2.304000000000000000000000000000000000000d-2


           x_star(2)=0.2304D-1


!        Column 3 of x_star
!
!         X_STAR(3)  = 6.D0/125.D0 * DSQRT(15.D0) - 47.D0/1875.D0
!
!          x_star(3)=0.1608365339512893d0

!          x_star(3) =1.608365339512893438219380725228885146533d-1


           x_star(3)=0.1608365339512893D0


!
!        Column 4 of x_star
!
!         X_STAR(4) = 488.D0/9375.D0

!          x_star(4)=0.5205333333333333d-1


           x_star(4)=0.5205333333333333D-1


!          x_star(4)=5.205333333333333333333333333333333333333d-2
!
!        Column 5 of x_star
!
!         X_STAR(5) = -6.D0/125.D0 * DSQRT(15.D0) - 47.D0/1875.D0
!

!          x_star(5)=-2.109698672846227d0d-1

!           x_star(5)=-0.2109698672846226771552714058562218479867


           x_star(5)=-0.2109698672846227D0


!        Column 6 of x_star
!
!         X_STAR(6) = 0.D0

          x_star(6)=0.0d0

!
       END IF
!
       IF (METHOD .EQ. 4) THEN
!
!        There are three new stages.
!
         S_STAR = 9
!
!        All these coefficients have been computed to a full 16 digits of
!        accuracy in Maple and then transferred here.

!         c_star(1) = 1.d0/2.d0 + sqrt(7.0d0)/14.0d0

          C_STAR(1) = 0.6889822365046136D0

          C_STAR(2) = 1.D0/5.D0 

          C_STAR(3) = 4.D0/5.D0 
!
!         v_star(1) = 1.d0/2.d0 + sqrt(7.0d0)/14.0d0

          V_STAR(1) = 0.6889822365046136D0

          V_STAR(2) = 1.D0/5.D0 

          V_STAR(3) = 4.D0/5.D0 
!
!        Column 1 of x_star
!
!      x_star(1) = 3.D0/392.D0*dsqrt(7.D0)-3.D0/392.D0
!      x_star(2) = -992.D0/140625.D0+152.D0/140625.D0*dsqrt(7.D0)
!      x_star(3) = 628.D0/140625.D0+152.D0/140625.D0*dsqrt(7.D0)
!
!        
!        Column 2 of x_star
!
!      x_star(4) = 3.D0/392.D0*dsqrt(7.D0)+3.D0/392.D0
!      x_star(5) = -628.D0/140625.D0+152.D0/140625.D0*dsqrt(7.D0)
!      x_star(6) = 992.D0/140625.D0+152.D0/140625.D0*dsqrt(7.D0)
!
!
!        Column 3 of x_star
!
!      x_star(7) = -5.D0/1176.D0*dsqrt(30.D0)*dsqrt(7.D0)
!     #-3.D0/10976.D0*dsqrt(30
!     #.D0)*dsqrt(525+70*dsqrt(30.D0))-3.D0/784.D0*dsqrt(7.D0)+27.D0/5488
!     #.D0*dsqrt(525+70*dsqrt(30.D0))
!      x_star(8) = 452.D0/109375.D0*dsqrt(525+70*dsqrt(30.D0))
!     #+193.D0/18750.D0*d
!     #sqrt(30.D0)+152.D0/140625.D0*dsqrt(30.D0)*dsqrt(7.D0)+9.D0/3125.D0
!     #-76.D0/9375.D0*dsqrt(7.D0)-226.D0/984375.D0*dsqrt(30.D0)*dsqrt(525
!     #+70*dsqrt(30.D0))
!      x_star(9) = -9.D0/3125.D0-76.D0/9375.D0*dsqrt(7.D0)
!     #+152.D0/140625.D0*dsqr
!     #t(30.D0)*dsqrt(7.D0)-193.D0/18750.D0*dsqrt(30.D0)-226.D0/984375.D0
!     #*dsqrt(30.D0)*dsqrt(525+70*dsqrt(30.D0))+452.D0/109375.D0*dsqrt(52
!     #5+70*dsqrt(30.D0))
!
!
!        Column 4 of x_star
!
!      x_star(10) = 5.D0/1176.D0*dsqrt(30.D0)*dsqrt(7.D0)
!     #+3.D0/10976.D0*dsqrt(525
!     #-70*dsqrt(30.D0))*dsqrt(30.D0)-3.D0/784.D0*dsqrt(7.D0)+27.D0/5488.
!     #D0*dsqrt(525-70*dsqrt(30.D0))
!      x_star(11) = 9.D0/3125.D0-76.D0/9375.D0*dsqrt(7.D0)
!     #-152.D0/140625.D0*dsqrt
!     #(30.D0)*dsqrt(7.D0)+226.D0/984375.D0*dsqrt(525-70*dsqrt(30.D0))*ds
!     #qrt(30.D0)-193.D0/18750.D0*dsqrt(30.D0)+452.D0/109375.D0*dsqrt(525
!     #-70*dsqrt(30.D0))
!      x_star(12) = -76.D0/9375.D0*dsqrt(7.D0)+452.D0/109375.D0
!     #*dsqrt(525-70*dsqr
!     #t(30.D0))+193.D0/18750.D0*dsqrt(30.D0)-9.D0/3125.D0-152.D0/140625.
!     #D0*dsqrt(30.D0)*dsqrt(7.D0)+226.D0/984375.D0*dsqrt(525-70*dsqrt(30
!     #.D0))*dsqrt(30.D0)
!
!
!        Column 5 of x_star
!
!      x_star(13) = 5.D0/1176.D0*dsqrt(30.D0)*dsqrt(7.D0)
!     #-3.D0/10976.D0*dsqrt(525
!     #-70*dsqrt(30.D0))*dsqrt(30.D0)-3.D0/784.D0*dsqrt(7.D0)-27.D0/5488.
!     #D0*dsqrt(525-70*dsqrt(30.D0))
!      x_star(14) = 9.D0/3125.D0-76.D0/9375.D0*dsqrt(7.D0)
!     #-152.D0/140625.D0*dsqrt
!     #(30.D0)*dsqrt(7.D0)-226.D0/984375.D0*dsqrt(525-70*dsqrt(30.D0))*ds
!     #qrt(30.D0)-193.D0/18750.D0*dsqrt(30.D0)-452.D0/109375.D0*dsqrt(525
!     #-70*dsqrt(30.D0))
!      x_star(15) = -76.D0/9375.D0*dsqrt(7.D0)-452.D0/109375.D0
!     #*dsqrt(525-70*dsqr
!     #t(30.D0))+193.D0/18750.D0*dsqrt(30.D0)-9.D0/3125.D0-152.D0/140625.
!     #D0*dsqrt(30.D0)*dsqrt(7.D0)-226.D0/984375.D0*dsqrt(525-70*dsqrt(30
!     #.D0))*dsqrt(30.D0)
!
!
!        Column 6 of x_star
!
!      x_star(16) = -5.D0/1176.D0*dsqrt(30.D0)*dsqrt(7.D0)
!     #+3.D0/10976.D0*dsqrt(30
!     #.D0)*dsqrt(525+70*dsqrt(30.D0))-3.D0/784.D0*dsqrt(7.D0)-27.D0/5488
!     #.D0*dsqrt(525+70*dsqrt(30.D0))
!      x_star(17) = 9.D0/3125.D0-76.D0/9375.D0*dsqrt(7.D0)
!     #+152.D0/140625.D0*dsqrt
!     #(30.D0)*dsqrt(7.D0)+193.D0/18750.D0*dsqrt(30.D0)+226.D0/984375.D0*
!     #dsqrt(30.D0)*dsqrt(525+70*dsqrt(30.D0))-452.D0/109375.D0*dsqrt(525
!     #+70*dsqrt(30.D0))
!      x_star(18) = -76.D0/9375.D0*dsqrt(7.D0)-193.D0/18750.D0
!     #*dsqrt(30.D0)-9.D0/
!     #3125.D0+152.D0/140625.D0*dsqrt(30.D0)*dsqrt(7.D0)+226.D0/984375.D0
!     #*dsqrt(30.D0)*dsqrt(525+70*dsqrt(30.D0))-452.D0/109375.D0*dsqrt(52
!     #5+70*dsqrt(30.D0))
!
!        Column 7 of x_star
!      x_star(19) = 0.0d0
!      x_star(20) = 4256.D0/140625.D0*dsqrt(7.D0)
!      x_star(21) = 4256.D0/140625.D0*dsqrt(7.D0)
!
!

        X_STAR(1) = 0.1259503554386166D-1
        X_STAR(2) = -0.4194459027329296D-2
        X_STAR(3) = 0.7325540972670704D-2
        X_STAR(4) = 0.2790115799284125D-1
        X_STAR(5) = -0.1606014582884851D-2
        X_STAR(6) = 0.9913985417115149D-2
        X_STAR(7) = 0.3142445823600003D-1
        X_STAR(8) = 0.140128203528661D0
        X_STAR(9) = 0.2161038635693079D-1
        X_STAR(10) = 0.1278455645496104D0
        X_STAR(11) = -0.2647240969774761D-1
        X_STAR(12) = 0.8052540747398258D-1
        X_STAR(13) = -0.2486766857825578D-1
        X_STAR(14) = -0.1547489916166944D0
        X_STAR(15) = -0.4775117444496425D-1
        X_STAR(16) = -0.1748985477440576D0
        X_STAR(17) = -0.3317969806100672D-1
        X_STAR(18) = -0.1516975152327369D0
        X_STAR(19) = 0.0D0
        X_STAR(20) = 0.8007336945700194D-1
        X_STAR(21) = 0.8007336945700194D-1

!        Column 8 of x_star
       X_STAR(22) = 0.0D0
       X_STAR(23) = 0.0D0
       X_STAR(24) = 0.0D0
!
!        Column 9 of x_star
       X_STAR(25) = 0.0D0
       X_STAR(26) = 0.0D0
       X_STAR(27) = 0.0D0

      END IF
!
      RETURN
      END SUBROUTINE TABLEAU_STAR


      SUBROUTINE COL_WEIGHTS(S_STAR,T,W)
!
! ***OVERVIEW***
!
!          This routine evaluates the weight polynomials (and their first 
!      derivatives) corresponding to the 's_star' stages used to form the 
!      interpolant for the discrete formula defined by the 'RK_tableau' 
!      routine. The evaluation takes place at 't'. The results are returned 
!      in the 's_star' locations of the vectors 'w' and 'wp'.  
!                            
!          The presence of the interp_weights routine demonstrates the use 
!      of an "information hiding" principle. The more common approach would be
!      to simply provide the coefficients of the polynomial weights in
!      an array made available to the rest of the software package through
!      the /interp_coeffs/ common block. In the approach employed here we
!      'hide' the polynomial expression for each weight inside the  
!      'interp_weights' routine whose function is to provide an evaluation
!      of the weights at a given point. An advantage of this approach is 
!      that changes can be made in the method by which the weights are 
!      determined with no subsequent change needed in the rest of the 
!      software. 
!
!----------------------------------------------------------------------------
!
!  ***declaration of parameters***
!
!      imports:
!
         INTEGER S_STAR
         DOUBLE PRECISION T
!
!        's_star' is the number of weight polynomials.
!        't'    is the point at which the weight polynomials 
!                 are to be evaluated. 
!
!      exports:
!
         DOUBLE PRECISION W(S_STAR)
!
!        'w'    contains the values of the weight polynomials
!               of the interpolant, evaluated at 't'.
!
!        'wp'   contains the values of the first derivative
!               of the weight polynomials of the interpolant,
!               evaluated  at 't'.
!




         INTEGER IMAXE

         COMMON /APPCOL/IMAXE




!------------------------------------------------------------------------------
!     called by: 'evaluate'
!------------------------------------------------------------------------------
!
!      The polynomials in 't' that define the interpolant are very
!      dependent on the underlying discrete formula and the  interpolant
!      derived. We use the value of 's_star' to determine which formula
!      is being implemented and thus which weights should be evaluated.
!
!      s_star = 4 implies number of collocation points is 2 and that
!      the CRK interpolant has four stages but b_4(t) = 0.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!  As part of SCSA technology, the more efficient approach is to hard code the weight polynomials for 1/2 and 2/3 across the subinterval
!
!  The only sampling pointsfor the error estimation in the ERRCHK routine are at 1/3 and 2/3 across the subinterval.
!



      IF (S_STAR .EQ. 4) THEN


        IF(IMAXE.NE.314)THEN



         IF(T.LT.0.2D0)THEN

             W(1)=0.9435739743742127D-1
           
             W(2)=-0.1764804844828675D-2

             W(3)=0.7407407407407407D-1

             W(4)=0.0D0


         ENDIF



         IF(T.LT.0.5D0.AND.T.GT.0.2D0)THEN

            W(1)=0.7330644658652999D-1

            W(2)=0.7676274875440806D-3

            W(3)=0.2592592592592593D0

            W(4)=0.0D0

          ENDIF


         IF(T.GT.0.5D0.AND.T.LT.0.8D0)THEN
     
            W(1)=-0.1397593589591207D0

            W(2)=0.6568528488504658D-1

            W(3)=0.7407407407407407D0

            W(4)=0.0D0


         ENDIF


         IF(T.GT.0.8D0)THEN

            W(1)=-0.2441436431433453D0

            W(2)=0.1515510505507527D0

            W(3)=0.9259259259259259D0

            W(4)=0.0D0

         ENDIF
    


        ENDIF


        IF(IMAXE.EQ.314)THEN

           W(1) = (3+DSQRT(3.D0))*(4*T**2-9*T+T*DSQRT(3.D0) &
                       +6-2*DSQRT(3.D0))*T/12


           W(2) = (3-DSQRT(3.D0))*(4*T-3+DSQRT(3.D0))*T**2/12


           W(3) = T**2*(3-2*T)


           W(4) = 0

        
         ENDIF




      END IF



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!      s_star = 6 implies number of collocation points is 3 and that
!      the CRK interpolant has six stages.

      IF (S_STAR .EQ. 6) THEN



        IF(IMAXE.NE.314)THEN



         IF(T.LT.0.2D0)THEN

!            w(1)= 5.679869684499314d-2

!             w(1)= 0.5679869684499314128943758573388203017833


            W(1)= 0.5679869684499314D-1


!            w(2)= -2.314814814814815d-2

!             w(2)= -0.2314814814814814814814814814814814814815

            W(2)=-0.2314814814814815D-1


!            w(3)= 9.491267653789124d-2

!            w(3)=  0.9491267653789125131376063850192043687221

            W(3)=0.9491267653789125D-1

!            w(4)= 1.186556927297668d-1

!            w(4)=0.1186556927297668038408779149519890260631


            W(4)=0.1186556927297668D0


!            w(5)= 5.340693937431719d-2

!            w(5)=0.5340693937431725348733675518806584570675

             W(5)=0.5340693937431725D-1
  

!            w(6)= -1.339591906721537d-1

!            w(6)= -0.1339591906721536351165980795610425240055


           W(6)=-0.1339591906721536D0


!          W(1) = T*(T-1)**2*(5*T**2-10*T+3)/3

!          W(2) = T**2*(T-1)*(30*T**2-35*T+9)/4

!          W(3) = (-5+DSQRT(15.D0))*(60*T**3-15*DSQRT(15.D0)*T**2-195* &
!               T**2+32*T*DSQRT(15.D0)+220*T-18*DSQRT(15.D0)-90)* &
!               T**2/36

!          W(4) = -4.D0/9.D0*T**2*(-18+62*T-75*T**2+30*T**3)

!          W(5) = -(5+DSQRT(15.D0))*(60*T**3-195*T**2+15*DSQRT(15.D0)* &
!                T**2+220*T-32*T*DSQRT(15.D0)+18*DSQRT(15.D0)-90)*  &
!                 T**2/36

!          W(6) = 125.D0/12.D0*T**2*(2*T-1)*(T-1)**2



         ENDIF




         IF(T.LT.0.5D0.AND.T.GT.0.2D0)THEN


!           W(1)=  1.097393689986283d-2

!            w(1)=  0.1097393689986282578875171467764060356653
          

             W(1)=0.1097393689986283D-1


!           W(2)=  -1.234567901234568

!            w(2)=-0.1234567901234567901234567901234567901235


           W(2)=-0.1234567901234568D-1


!           W(3)=  2.468582536682419

!            w(3)=0.2468582536682418898540264855907379865797

           W(3)=0.2468582536682419D0


!           W(4)=  2.249657064471879d-1

!            w(4)=0.2249657064471879286694101508916323731138


            W(4)=0.2249657064471879D0


!           W(5)=  3.434887939074224

!            w(5)=0.3434887939074302098273620302380247981295


            W(5)=0.3434887939074302D-1


!           W(6)=  -1.714677640603567


!            w(6)=-0.1714677640603566529492455418381344307271

            W(6)=-0.1714677640603567D0


!         W(1) = T*(T-1)**2*(5*T**2-10*T+3)/3


!         W(2) = T**2*(T-1)*(30*T**2-35*T+9)/4

!         W(3) = (-5+DSQRT(15.D0))*(60*T**3-15*DSQRT(15.D0)*T**2-195* &
!              T**2+32*T*DSQRT(15.D0)+220*T-18*DSQRT(15.D0)-90)* &
!              T**2/36

!         W(4) = -4.D0/9.D0*T**2*(-18+62*T-75*T**2+30*T**3)

!         W(5) = -(5+DSQRT(15.D0))*(60*T**3-195*T**2+15*DSQRT(15.D0)* &
!               T**2+220*T-32*T*DSQRT(15.D0)+18*DSQRT(15.D0)-90)*  &
!                T**2/36

!         W(6) = 125.D0/12.D0*T**2*(2*T-1)*(T-1)**2


        ENDIF




         IF(T.GT.0.5D0.AND.T.LT.0.8D0)THEN
     



!                W(1)=  -3.566529492455418d-2

!                 w(1)=-0.3566529492455418381344307270233196159122

                 W(1)=-0.3566529492455418D-1


!                W(2)=  3.703703703703704d-2

!                 w(2)=0.3703703703703703703703703703703703703704

                 W(2)=0.3703703703703704D-1



!               W(3)=  3.496835855257845d-1


                

!                 w(3)=0.3496835855257841912306867160374430513483


                 W(3)=0.3496835855257842D0


!               W(4)=  2.194787379972565d-1

!                 w(4)=0.2194787379972565157750342935528120713306


                 W(4)= 0.2194787379972565D0


!               W(5)=  -7.533516302921969d-2

!                w(5)=-0.7533516302921354651189384909642796218260

                 W(5)= -0.7533516302921355D-1

!               W(6)=  1.714677640603567d-1

!                w(6)=0.1714677640603566529492455418381344307271

                W(6)=0.1714677640603567D0



!         W(1) = T*(T-1)**2*(5*T**2-10*T+3)/3

!         W(2) = T**2*(T-1)*(30*T**2-35*T+9)/4

!         W(3) = (-5+DSQRT(15.D0))*(60*T**3-15*DSQRT(15.D0)*T**2-195* &
!              T**2+32*T*DSQRT(15.D0)+220*T-18*DSQRT(15.D0)-90)* &
!              T**2/36

!         W(4) = -4.D0/9.D0*T**2*(-18+62*T-75*T**2+30*T**3)

!         W(5) = -(5+DSQRT(15.D0))*(60*T**3-195*T**2+15*DSQRT(15.D0)* &
!               T**2+220*T-32*T*DSQRT(15.D0)+18*DSQRT(15.D0)-90)*  &
!                T**2/36

!         W(6) = 125.D0/12.D0*T**2*(2*T-1)*(T-1)**2







         ENDIF




         IF(T.GT.0.8D0)THEN

!          w(1)=-1.436042524005487

        


!           w(1)= -0.1436042524005486968449931412894375857339

            W(1)=-0.1436042524005487D-1


!          w(2)=-1.929012345679012

!           w(2)=-0.1929012345679012345679012345679012345679

            W(2)=-0.1929012345679012D-1


!          w(3)=3.073823127306108d-1

!           w(3)=0.3073823127306085199432887892174211144027

            W(3)=0.3073823127306085D0


!          w(4)=3.257887517146776

!           w(4)=0.3257887517146776406035665294924554183813

            W(4)=0.3257887517146776D0


!          w(5)=9.985362691272629

!           w(5)=0.9985362691273853081116937264814815858089


            W(5)=0.9985362691273853D-1


!          w(6)=1.339591906721537d-1

!           w(6)=0.1339591906721536351165980795610425240055


             W(6)=0.1339591906721536D0


!          W(1) = T*(T-1)**2*(5*T**2-10*T+3)/3

!          W(2) = T**2*(T-1)*(30*T**2-35*T+9)/4

!          W(3) = (-5+DSQRT(15.D0))*(60*T**3-15*DSQRT(15.D0)*T**2-195* &
!               T**2+32*T*DSQRT(15.D0)+220*T-18*DSQRT(15.D0)-90)* &
!               T**2/36

!          W(4) = -4.D0/9.D0*T**2*(-18+62*T-75*T**2+30*T**3)

!          W(5) = -(5+DSQRT(15.D0))*(60*T**3-195*T**2+15*DSQRT(15.D0)* &
!                T**2+220*T-32*T*DSQRT(15.D0)+18*DSQRT(15.D0)-90)*  &
!                 T**2/36

!          W(6) = 125.D0/12.D0*T**2*(2*T-1)*(T-1)**2




         ENDIF
    
        ENDIF

        IF(IMAXE.EQ.314)THEN

          
         W(1) = T*(T-1)**2*(5*T**2-10*T+3)/3

         W(2) = T**2*(T-1)*(30*T**2-35*T+9)/4

         W(3) = (-5+DSQRT(15.D0))*(60*T**3-15*DSQRT(15.D0)*T**2-195* &
              T**2+32*T*DSQRT(15.D0)+220*T-18*DSQRT(15.D0)-90)* &
              T**2/36

         W(4) = -4.D0/9.D0*T**2*(-18+62*T-75*T**2+30*T**3)

         W(5) = -(5+DSQRT(15.D0))*(60*T**3-195*T**2+15*DSQRT(15.D0)* &
               T**2+220*T-32*T*DSQRT(15.D0)+18*DSQRT(15.D0)-90)*  &
                T**2/36

         W(6) = 125.D0/12.D0*T**2*(2*T-1)*(T-1)**2
        

         ENDIF




!         W(1) = T*(T-1)**2*(5*T**2-10*T+3)/3

!         W(2) = T**2*(T-1)*(30*T**2-35*T+9)/4

!         W(3) = (-5+DSQRT(15.D0))*(60*T**3-15*DSQRT(15.D0)*T**2-195* &
!              T**2+32*T*DSQRT(15.D0)+220*T-18*DSQRT(15.D0)-90)* &
!              T**2/36

!         W(4) = -4.D0/9.D0*T**2*(-18+62*T-75*T**2+30*T**3)

!         W(5) = -(5+DSQRT(15.D0))*(60*T**3-195*T**2+15*DSQRT(15.D0)* &
!               T**2+220*T-32*T*DSQRT(15.D0)+18*DSQRT(15.D0)-90)*  &
!                T**2/36

!         W(6) = 125.D0/12.D0*T**2*(2*T-1)*(T-1)**2


      END IF


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


      IF (S_STAR .EQ. 9) THEN



        IF(IMAXE.NE.314)THEN



           IF(T.LT.0.2D0)THEN

              W(1)=-0.1150348945591317D0

              W(2)=0.1387650209689777D0

              W(3)= 0.4819275638059541D0

              W(4)=0.4118225626662212D0

              W(5)=-0.2296787116499782D0

              W(6)=-0.3847718606383842D0

              W(7)=-0.2281661873591156D0

              W(8)=-0.5024271855143026D0

              W(9)=0.5942303589464258D0



           ENDIF




           IF(T.LT.0.5D0.AND.T.GT.0.2D0)THEN

              W(1)= -0.4573080643257935D-1

              W(2)=0.7095561972323692D-1

              W(3)=0.2894461668828698D0

              W(4)=0.266482498906871D0

              W(5)=-0.9382952923818732D-1

              W(6)=-0.1973529088423723D0

              W(7)=-0.1446506068555011D0

              W(8)= -0.1205466853346099D0

              W(9)=0.3085595845236023D0





           ENDIF



           IF(T.GT.0.5D0.AND.T.LT.0.8D0)THEN
     
              W(1)=0.4961740908461741D-1

              W(2)=-0.7484222237526824D-1

              W(3)=0.3788614407585321D-1

              W(4)=0.1731351307224674D0

              W(5)=0.3063570544713461D0

              W(6)=0.2178754430209323D0

              W(7)=0.1446506068554992D0

              W(8)=0.1999299094367667D0

              W(9)=-0.3879428086256833D0




           ENDIF




            IF(T.GT.0.8D0)THEN

              W(1)=0.5142214021528287D-1

              W(2)=-0.7515226662512127D-1

              W(3)=0.3281638326198989D-1

              W(4)=0.1665107508495575D0

              W(5)=0.3034905529965711D0

              W(6)=0.21788275870745D0

              W(7)=0.2281661873590951D0

              W(8)=0.2078410023235988D0

              W(9)=-0.2996441757555419D0



            ENDIF
    


        ENDIF



        IF(IMAXE.EQ.314)THEN

      W(1)=1.D0*T-0.3059376195485321D2*T**2+0.1903161507410791D3*T**3 &
      -0.5066326577159306D3*T**4+0.6779074502047026D3*T**5- &
      0.4502127204516637D3*T**6+0.1182155391766658D3*T**7 

      W(2)=0.1826042862151987D2*T**2-0.1223369840744124D3*T**3 &
      +0.3167368243825973D3*T**4-0.401428283538036D3*T**5 &
      +0.2496918871183304D3*T**6-0.6092387250999916D2*T**7

      W(3)=0.6644072302985358D2*T**2-0.4659005706172424D3*T**3 &
      +0.1290052895837388D4*T**4-0.1760013679899507D4*T**5 &
      +0.1182484831564531D4*T**6-0.3128902724924552D3*T**7

      W(4)=0.5805733615326762D2*T**2-0.4153918067568109D3*T**3 &
      +0.1181297756609702D4*T**4-0.1652285224203471D4*T**5 &
      +0.113424023288076D4*T**6-0.3055922221060169D3*T**7

      W(5)=-0.2871080418445304D2*T**2+0.1822499138934506D3*T**3 &
      -0.4329611914049302D3*T**4+0.4970100823644704D3*T**5 &
      -0.2782997171236686D3*T**6+0.6103778903256205D2*T**7

      W(6)=-0.5078725499866815D2*T**2+0.3415424634806026D3*T**3 &
      -0.8908894610421601D3*T**4+0.1143788821738507D4*T**5 &
      -0.7259253473216229D3*T**6+0.18244470556591D3*T**7

      W(7)=-0.3266666666666667D2*T**2+0.2381228070175439D3*T**3 &
      -0.7006140350877193D3*T**4+0.1032438596491228D4*T**5 &
      -0.7521929824561404D3*T**6+0.2149122807017544D3*T**7

      W(8)=-0.7911786699358921D2*T**2+0.5889108452680835D3*T**3 &
      -0.1685448258473616D4*T**4+0.2349773462964052D4*T**5 &
      -0.1603256196611642D4*T**6+0.4291380138467126D3*T**7

      W(9)=0.7911786699358921D2*T**2-0.537512818952294D3*T**3 &
      +0.1428458126894669D4*T**4-0.1887191226121947D4*T**5 &
      +0.1243470012401116D4*T**6-0.3263419612151337D3*T**7



         ENDIF



      END IF



      RETURN
      END SUBROUTINE COL_WEIGHTS



       SUBROUTINE FSUB_FIRST(NODE,NSOL,ORDERS,Z,F,F_FIRST)
        !The purpose of this routine is to compute the function evaluations for a first
        !order set of differential equations given only a system of mixed order equations.
        !This requires as input the number of differential equations of the mixed order
        ! system, NODE. It also requires NSOL solution and derivative components stored in Z,
        ! the order of each differential equation in ORDERS and the function evaluations, F,
        ! from the previous call to fsub (mixed order routine). The results are stored in the
        ! NSOL components of F_FIRST.
        !------------------------------------------------------------------------------------
        ! CALLED BY: SETUP_PHASE1, SETUP
        !----------------------------------------------------------------------------------- 
        ! Input Arguments
        INTEGER  NODE,NSOL
        DOUBLE PRECISION Z(NSOL)
        DOUBLE PRECISION F(NODE)
        INTEGER ORDERS(NODE)
        ! Output Arguments
        DOUBLE PRECISION F_FIRST(NSOL)
        ! Local Variables
        INTEGER I,J 
        INTEGER ORDERS_SUM 
       
        ! Store the function evaluations from the previous call to fsub in the appropriate
        ! locations of F_FIRST. Then the remaining equations are stored as F_FIRST(I)=Z(I+1).
        ! The remaining equations are of the form y_i'=y_i+1.


        F_FIRST(ORDERS(1))=F(1)

        DO I=1,ORDERS(1)-1
          F_FIRST(I)=Z(I+1)
        ENDDO

        ORDERS_SUM=ORDERS(1)

        DO J=2,SIZE(ORDERS)
          ORDERS_SUM=ORDERS_SUM+ORDERS(J)
          F_FIRST(ORDERS_SUM)=F(J)

          DO I=(ORDERS_SUM-ORDERS(J)+1),ORDERS_SUM-1
             F_FIRST(I)=Z(I+1)
          ENDDO

        ENDDO

       RETURN
       END SUBROUTINE FSUB_FIRST


       SUBROUTINE SETUP_PHASE1(NSUB,MESH,ASAVE,ZVAL,DMZ,K_POINTS,NCOMP,&
                  MMAX,M,MSTAR,Y_SUPER,SK_DISCRETE_SUPER, &
                  SK_INTERP_SUPER,FSUB)                

        IMPLICIT NONE

        ! See page 233 of Enright and Muir 1999, Superconvergent Interpolants for the 
        ! Collocation Solution of BVODES where it states the purpose of this routine.
        ! It states a superinterpolant is obtained from a CRK scheme which, on the 
        ! ith subinterval will use the superconvergent endpoints values y_i-1 and y_i,
        ! the corresponding stages f(t_i-1,y_i-1) and f(t_i,y_i), the collocation stages,
        ! and as many extra MIRK stages as necessary to achieve the desired order.
     
        ! Parameter values to change the size of the arrays to allow for more than
        ! ( mstar = m(1) + ... + m(ncomp) .le. 40 )
        INTEGER MAXMSTAR
        PARAMETER (MAXMSTAR=160)
        INTEGER MAXNODE
        PARAMETER (MAXNODE=50)

        ! Input Variables
        INTEGER K_POINTS ! k collocation points
        INTEGER NSUB ! number of subintervals
        DOUBLE PRECISION MESH(NSUB+1) ! mesh points
        DOUBLE PRECISION ASAVE(7,4) ! Dummy variable-maybe removed later
        DOUBLE PRECISION ZVAL(MAXMSTAR*(NSUB+1)) ! The current solution approximation at the mesh points. 
        DOUBLE PRECISION DMZ(K_POINTS*MAXNODE*NSUB) ! To be described later
        INTEGER NCOMP ! The number of differential equations
        INTEGER MMAX ! The highest order of the differential equations
        INTEGER M(NCOMP) ! The order of each differential equation
        INTEGER MSTAR ! The sum of the orders in the differential eqation, also the number of solution components of the Z vector

        ! Output Variables
        DOUBLE PRECISION Y_SUPER(MAXMSTAR*(NSUB+1))
        DOUBLE PRECISION SK_DISCRETE_SUPER(NSUB*10*MAXMSTAR)
        DOUBLE PRECISION SK_INTERP_SUPER(NSUB*10*MAXMSTAR)
        
        ! Local Variables
        INTEGER S ! Number of collocation stages plus the two stages at the ends of a subinterval
        INTEGER S_STAR ! Number of collocation stages plus the two stages at the ends of a subinterval plus the extra MIRK stages. 
        DOUBLE PRECISION XC(4) ! The C values in the tableau for the collocation stages
        INTEGER I ! Loop counter
        DOUBLE PRECISION X ! Represents the mesh point
        DOUBLE PRECISION Z(MAXMSTAR) ! Approximate solution at the mesh point
        DOUBLE PRECISION DUMMY(1) ! Dummy variable - maybe deleted later
        INTEGER K ! Loop counter
        DOUBLE PRECISION F(MAXNODE) ! The mixed order fsub evaluations
        DOUBLE PRECISION F_FIRST(MAXMSTAR) ! The first order equivalent of the fsub evaluations
        DOUBLE PRECISION HI ! Length of the subinterval
        INTEGER L ! Loop counter - maybe deleted later
        INTEGER ITEST,JTEST ! Integer needed for call to approx
        INTEGER J,IZ ! Loop counter
        DOUBLE PRECISION YY(MAXMSTAR) ! Local variable needed for call to setup
        DOUBLE PRECISION A(28) ! Local varible needed for call to approx

        ! Input constants for approx routine through common block
        DOUBLE PRECISION RHO(7),COEF(49)
        COMMON /COLLOC/ RHO,COEF

        EXTERNAL FSUB


        ! The first 2 stages are f(t_i-1,y_i-1) and f(t_i,y_i) and next k stages are the
        ! collocation stages hence we have the number k_points+2 for S.
        ! The extra MIRK stages are 0 for k=2, 1 for k=3, and 3 for k=4.
        S = K_POINTS+2
        IF (K_POINTS .EQ. 2) S_STAR = K_POINTS+2
        IF (K_POINTS .EQ. 3) S_STAR = K_POINTS+2+1
        IF (K_POINTS .EQ. 4)  S_STAR =  K_POINTS+2+3

        !The main purpose of this routine is to define the parameter which
        !gives,'s_star', the total number of stages for the interpolant and the 
        !coefficients for the 's_star-s' extra stages needed by the interpolant.
        !There is a common block inside the routine tableau_star which passes these
        !coefficients for c,v and a for the extra MIRK stages.

        CALL TABLEAU_STAR(K_POINTS)


        ! Following are the C values in the tableaus for the 1999 Enright and muir paper,
        ! corresponding to the collocation stages. Currently only tableaus for k=2,3 and 4
        ! have been determined. They are called XC here because they are special C values in
        ! in the tableau for the collocation stages. They represent the relative distance across
        ! each subinterval. 

        ! Hard code values for XC to minimize cost of SCSA technology

        IF (K_POINTS .EQ. 2) THEN

!               XC(1) = 1.0D0/2.0D0 - SQRT(3.0D0)/6.0D0

                XC(1)= 0.2113248654051871D0

!               XC(2) = 1.0D0/2.0D0 + SQRT(3.0D0)/6.0D0
                  
                XC(2)= 0.7886751345948129D0


        END IF

        IF (K_POINTS .EQ. 3) THEN


!                XC(1) = 1.0D0/2.0D0 - SQRT(15.0D0)/10.0D0

                 XC(1)= 0.1127016653792583D0

!                XC(2)=1.0D0/2.0D0

                 XC(2) = 0.5d0 

!                XC(3) = 1.0D0/2.0D0 + SQRT(15.0D0)/10.0D0

                 XC(3)= 0.8872983346207417D0




        END IF

        IF (K_POINTS .EQ. 4) THEN

!               XC(1) = 1.0D0/2.0D0 - SQRT(525.0D0+70.0D0*SQRT(30.0D0))/70.0D0

                XC(1)= 0.6943184420297371D-1

!               XC(2) = 1.0D0/2.0D0 - SQRT(525.0D0-70.0D0*SQRT(30.0D0))/70.0D0

                XC(2)=0.3300094782075719D0

!               XC(3) = 1.0D0/2.0D0 + SQRT(525.0D0-70.0D0*SQRT(30.0D0))/70.0D0

                XC(3)= 0.6699905217924281D0

!               XC(4) = 1.0D0/2.0D0 + SQRT(525.0D0+70.0D0*SQRT(30.0D0))/70.0D0

                XC(4)= 0.9305681557970263D0


        END IF


        ! For each subinterval compute corresponding stages f(t_i-1,y_i-1) and  
        ! f(t_i,y_i).

        DO  I=1,(NSUB+1)

          ! x=t_i
          X=MESH(I)


           
            DO  J = 1, MSTAR
              Z(J) = ZVAL((I-1)*MSTAR+J)
              Y_SUPER((I-1)*MSTAR+J)=Z(J)
            ENDDO


          ! Compute the f(t_i,y_i) stage with the call to fsub mixed order
          CALL FSUB(X,Z,F)

          ! In the 1999 Enright and Muir paper the function evaluations have to be
          ! in first order form. f_first is then returned with nsol or mstar solution components. 
          

          ! Store the stages in sk_discete_super  
          IF (I .NE. 1) THEN
            DO K = 1, MSTAR
              SK_DISCRETE_SUPER((I-2)*S*MSTAR + MSTAR + K)=F(K)
            ENDDO
          END IF

          IF (I .NE. (NSUB+1)) THEN
            DO K = 1, MSTAR
              SK_DISCRETE_SUPER((I-1)*S*MSTAR + K)=F(K)
            ENDDO
          END IF

        ENDDO ! End of do i=1,Nsub+1 loop

             
        ! Now compute the k collocation stages on each subinterval

        DO  I= 1, NSUB

          ! Compute hi, the distance between mesh points, and use xc to
          ! find the distance across each subinterval aand store in x.
          HI = MESH(I+1)-MESH(I)
          DO K=1,K_POINTS

            X = MESH(I)+HI*XC(K)
            ! Compute y(t_i+c*hi) with the cal to approx stored in z.

            CALL APPROX(I,X,Z,A,COEF,MESH,NSUB,ZVAL,DMZ,K_POINTS,&
                 NCOMP,MMAX,M,MSTAR,3,DUMMY,0)
  
            ! Compute the f(t_i+c*hi,y(t_i+c*hi) stage with the call to fsub mixed order
            CALL FSUB(X,Z,F)
 
            ! In the 1999 Enright and Muir paper the function evaluations have to be
            ! in first order form. f_first is then returned with nsol or mstar solution components.

            ! Store the collocation stages in sk_discrete_super.

            DO L = 1, MSTAR
              SK_DISCRETE_SUPER((I-1)*S*   &
              MSTAR+(2+K-1)*MSTAR+L)=F(L)
            ENDDO

           ENDDO ! End of do k=1,k_points loop

         ENDDO ! End of do i=1,Nsub loop


         ! Compute the extra MIRK stages for k=3 and k=4.
         IF (K_POINTS .NE. 2) THEN

           ! Compute the extra MIRK stages on each subinterval.
           ! Store them in sk_interp_super.

           DO I=1,NSUB
     
           CALL SETUP(NCOMP,MSTAR,M, &
           MESH(I),MESH(I+1)-MESH(I), &
           ZVAL((I-1)*MSTAR+1),ZVAL(I*MSTAR+1), S, &
           SK_DISCRETE_SUPER((I-1)*MSTAR*(K_POINTS+2)+1), S_STAR, &
           SK_INTERP_SUPER((I-1)*MSTAR*(S_STAR-K_POINTS-2)+1),YY,FSUB)

           ENDDO
         END IF


         RETURN
         END SUBROUTINE SETUP_PHASE1


END MODULE COLNEW_M
