subroutine movcol(npde, npts, atol, rtol, touta, ntouta, par,
     +                  iflag, rwork, lrw, iwork, liw)
c
c version: Nov. 22, 1995 
c
c 1/17/96  change integer phypde to logical phypde
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                     c
c PURPOSE :                                                           c
c                                                                     c
c          MOVCOL is primarily intended to solve systems of           c
c          second-order parabolic PDEs in one space dimension. It     c
c          is also capable of solving hyperbolic PDEs with suitably   c
c          smooth solutions [HR1].                                    c
c                                                                     c
c          MOVCOL uses a method of lines approach based upon a        c
c          MOVing COLlocation method [HR1].  The physical PDEs are    c
c          discretized in space with a cubic Hermite colloction-      c
c          type method, and the MMPDEs (moving mesh PDEs) for         c
c          computing the moving mesh points [HRR1] are discretized    c
c          with a 3-point finite difference method.  The resulting    c
c          ODE system is integrated in time with the DAE solver       c
c          DASSL developed by L. Petzold [Petzold].                   c 
c                                                                     c
c          The physical PDEs are considered in the general divergence c
c          form                                                       c
c                                                                     c
c          F(t, x, u, ux, ut, uxt) = [ G(t, x, u, ux, ut, uxt) ]_x    c
c                                                                     c
c                                            x^L(t) < x < x^R(t)      c
c                                            t_a < t < t_b            c
c                                                                     c
c          subject to the separated boundary conditions               c
c                                                                     c
c          B^L(t, x, xt, u, ux, uxx, ut, uxt) = 0   at x = x^L(t)     c
c          B^R(t, x, xt, u, ux, uxx, ut, uxt) = 0   at x = x^R(t)     c
c                                                                     c
c          and the initial conditions                                 c
c                                                                     c
c          u(x,t_a) = U(x)                                            c
c                                                                     c
c          where F, G, B^L, B^R, u and U are vector-valued functions  c
c          of length npde. It is assumed that the boundaries x^L(t)   c
c          and x^R(t) are determined by the two conditions            c
c                                                                     c
c          b^L_M(t, x, xt, u, ux, uxx, ut, uxt) = 0   at x = x^L(t)   c
c          b^R_M(t, x, xt, u, ux, uxx, ut, uxt) = 0   at x = x^R(t)   c
c                                                                     c
c          where b^L_M and b^R_M are scalar functions.                c
c                                                                     c
c          It is the user's responsibility to ensure that the PDEs    c
c          are well posed.                                            c
c                                                                     c
c          Notation :                                                 c
c                                                                     c
c          t      -- time                                             c
c          x      -- physical spatial coordinate                      c
c          x^L(t) -- left boundary                                    c
c          x^R(t) -- right boundary                                   c
c          u      -- physical solution variable                       c
c          xt     -- mesh speed x_t                                   c
c          ux     -- 1st derivative of u with respect to x            c
c          uxx    -- 2nd derivative of u with respect to x            c
c          ut     -- 1st derivative of u with respect to t            c
c          uxt    -- 2nd mixed derivative                             c
c                                                                     c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                     c
c AUTHORS:                                                            c
c                                                                     c
c          Weizhang Huang                                             c
c          Department of Mathematics                                  c
c          the University of Kansas                                   c
c          Lawrence, KS 66044                                         c
c          USA                                                        c
c          whuang@math.ku.edu                                      c
c                                                                     c
c          Robert D. Russell                                          c
c          Department of Mathematics and Statistics                   c
c          Simon Fraser University                                    c
c          Burnaby, B.C. V5A 1S6, Canada                              c
c          rdr@cs.sfu.ca                                              c
c                                                                     c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                     c
c REFERENCES :                                                        c
c                                                                     c
c 1. [HR1] W. Huang and R. D. Russell,                                c
c    A moving collocation method for solving time dependent partial   c
c    differential equations,                                          c
c    Appl. Numer. Math. 20 (1996), 101-116.                           c
c 2. [HR2] W. Huang and R. D. Russell,                                c
c    Analysis of moving mesh partial differential equations with      c
c    spatial smoothing,                                               c
c    SIAM J. Numer. Anal. 34 (1997), 1106-1126.                       c
c 3. [HRR1] W. Huang, Y. Ren and R. D. Russell,                       c
c    Moving mesh partial differential equations (MMPDEs) based on the c
c    equidistribution principle,                                      c
c    SIAM J. Numer. Anal. 31 (1994), 709--730.                        c
c 4. [HRR2] W. Huang, Y. Ren and R. D. Russell,                       c
c    Moving mesh methods based on moving mesh partial differential    c
c    equations,                                                       c
c    J. Comput. Phys. 113 (1994), 279--290.                           c
c 5. [Petzold] L. R. Petzold,                                         c
c    A description of DASSL: A differential/algebraic system solver,  c
c    SAND 82-8637, Sandia Labs, Livermore, USA                        c
c                                                                     c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                     c
c USAGE :                                                             c
c                                                                     c
c      integer          npde, npts, ntouta, iflag, lrw, liw,          c
c     +                 iwork(liw)                                    c
c      double precision atol, rtol, touta(ntouta), par(20),           c
c     +                 rwork(lrw)                                    c
c                                                                     c
c      call movcol(npde, npts, atol, rtol, touta, ntouta, par,        c
c     +            iflag, rwork, lrw, iwork, liw)                     c
c                                                                     c
c      NOTE :                                                         c
c                                                                     c
c      lrw  >= 60 + 6*npde + (22+3*ip)*npts                           c
c            + (55+12*ip)*npts*npde + (24+12*ip)*npts*npde**2         c
c      liw  >= 20 + npts + 2*npts*npde                                c
c                                                                     c
c      where ip := par(6)                                             c
c                                                                     c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                     c
c DESCRIPTION OF VARIABLES ( I : input, O : output)                   c
c                                                                     c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c 1. description of parameters
c
c npde:I    (integer) number of physical pdes
c
c npts:I    (integer) number of mesh nodes
c
c atol:I    (real) absolute tolerance for time integration 
c
c rtol:I    (real) relative tolerance for time integration 
c
c touta:I   (real array) prescribed times for outputting information
c           about performance of time integration with DASSL.
c           the ode system is integrated from touta(1) := t_a to
c           touta(ntouta) := t_b
c
c ntouta:I  (integer) dimension of touta
c
c par:I     (real array) array of length 20 containing optional
c           parameters
c
c   NOTE:   if a parameter does not have a value within its permitted
c           range, then its default value is used. therefore, if a
c           default value is desired, the parameter can be set to a
c           value outside of its permitted range, for example, to a 
c           negative number (such as -1) 
c
c   par  name     meaning
c
c   1   job   - parameter specifying whether or not an equidistributed
c                mesh is to be generated with respect to the initial 
c                solution
c
c                job = 1: solve the physical pdes starting with a mesh
c                         provided by the user in the subroutine DEFMSH
c                job = 2: solve the physical pdes after first computing
c                         an equidistributed mesh
c                job = 3: generate an equidistributed mesh and stop
c
c                NOTE: in the cases job = 2 and 3, the physical region
c                must be defined by setting par(8) = x^L(t_a) and
c                par(9) = x^R(t_a). the mesh generation uses the
c                initial value U(x) of the physical solution provided
c                by the user in the subroutine DEFIVS. the mesh is
c                generated internally by solving the user-prescribed
c                MMPDE and the pde u_t = U(x), 0 <= t <= 1, with the
c                user-provided monitor function
c
c                PERMITTED VALUES: 1 <= job <= 3
c
c                DEFAULT VALUE: job = 2 
c
c   2   inform - parameter for selection of output information about
c                performance of the time integration code DASSL 
c
c                PERMITTED VALUES: 1 <= inform 
c
c                DEFAULT VALUE: inform = 6
c 
c   3   mmpde  - parameter used for specifying the basic moving mesh
c                pde 
c
c                mmpde = 1: fixed mesh
c                mmpde = 2: MMPDE4 in [HRR1, HRR2]
c                mmpde = 3: MMPDE6 in [HRR1, HRR2]
c                mmpde = 4: MMPDE(24) with spatial smoothing in [HR2],
c                           a smoothed version of MMPDE4 in terms of
c                           the mesh concentration function
c
c                PERMITTED VALUES: 1 <= mmpde <= 4
c
c                DEFAULT VALUE: mmpde = 4 
c
c   4   tau    - time smoothing parameter in MMPDEs
c
c                PERMITTED VALUES: 0 <= tau
c
c                RECOMMENDED VALUES: 1e-6 <= tau <= 1e-3
c
c                DEFAULT VALUE: tau = 1e-4
c
c   5   gamma  - spatial smoothing parameter used when mmpde = 4
c
c                PERMITTED VALUES: 0 < gamma
c
c                RECOMMENDED VALUES: 1 <= gamma <= 2
c
c                DEFAULT VALUE: gamma = 1 
c
c   6   ip     - spatial smoothing parameter used when mmpde = 2 or 3
c
c                PERMITTED VALUES: 0 <= ip
c 
c                RECOMMENDED VALUES: 2 <= ip <= 4
c
c                DEFAULT VALUE: ip = 2
c
c   7   stpmax - the allowed maximal time stepsize for the time
c                integration
c
c                PERMITTED VALUES: 0 < stpmax 
c
c                DEFAULT VALUE: no restriction on time stepsize
c
c   8   xl     - the left end of the initial phyical region, i.e. xl =
c                x_L(t_a). this value is used when job = 2 and 3
c
c   9   xr     - the right end of the initial phyical region, i.e. xr =
c                x_R(t_a). this value is used when job = 2 and 3
c
c   10 - 20    - not used 
c
c iflag:O   (integer) an indicator reporting the code's performance 
c
c           iflag = +1:  successful run and default values of some of
c                        parameters par(1) -- par(7) have been used
c           iflag =  0:  successful run
c           iflag = -1:  at least one of the workspace arrays, rwork
c                        or iwork, has a length less than that required
c           iflag = -2:  at least one of the parameters, npde, npts,
c                        atol, rtol or ntouta, is incorrectly input
c           iflag = -3:  time integration stops since istop is set
c                        less than zero by the user in DEFOUT 
c           iflag = -4:  time integration stops since IDID <= -2 in
c                        DASSL (see documentation in DASSL) 
c
c 2. workspace arrays
c
c rwork     (real array) workspace 
c
c lrw:I     (integer) dimension of rwork. it must be at least
c           60 + 6*npde + (22+3*ip)*npts + (55+12*ip)*npts*npde
c           + (24+12*ip)*npts*npde**2, where ip := par(6)
c           (if ipar(3)=mmpde.ne.2 or 3, then one can consider
c           ip to be 0 in this lower bound for lrw)
c
c iwork     (integer array) workspace 
c
c liw:I     (integer) dimension of iwork. it must be at least
c           20 + npts + 2*npts*npde
c
c 3. subroutines for problem definition
c
c DEFPDE    define the physical pdes
c DEFBCP    define the physical boundary conditions
c DEFIVS    define the initial values for the physical solution and
c           its first spatial derivative 
c DEFMSH    define the initial mesh
c DEFBCM    define the boundary conditions for the mesh 
c DEFMNT    define the monitor function
c
c see section 5 below for a detailed description of these subroutines
c
c 4. subroutines for solution output 
c
c DEFOUT    output the physical solution and mesh
c DEFFNC    interpolate the physical solutions and first spatial
c           derivative via cubic Hermite interpolation. this is
c           provided for user evaluation of solutions at non-mesh
c           points
c
c see section 5 below for a detailed description of these subroutines
c
c 5. subroutines (I: included, U: user must provide)
c
c    NOTE : no input variable values should be changed in the user
c           provided subroutines 
c
c DEFPDE:U  subroutine DEFPDE (npde, index, t, x, u, ux, ut, uxt, fg)
c           integer          npde, index
c           double precision t, x, u(npde), ux(npde), ut(npde),
c                            uxt(npde), fg(npde)
c
c           this subroutine is used to define the physical pdes. the F
c           and G terms, given by fg, must be defined in terms of input
c           variables npde, t, x, u, ux, ut and uxt by
c
c           fg := F(t, x, u, ux, ut, uxt)      if index < 0 
c           fg := G(t, x, u, ux, ut, uxt)      if index > 0 
c
c DEFBCP:U  subroutine DEFBCP(npde, index, t, x, xt, u, ux, uxx, ut,
c                             uxt, res)
c           integer          npde, index
c           double precision t, x, xt, u(npde), ux(npde), uxx(npde),
c                            ut(npde), uxt(npde), res(npde)
c
c           this subroutine is used to define the boundary conditions
c           for the physical solution variables. the residual res must
c           be defined in terms of input variables npde, index, t, x,
c           xt, u, ux uxx, ut and uxt by
c
c           res(k) := B^L_k(t, x, xt, u, ux, uxx, ut, uxt) if index < 0
c           res(k) := B^R_k(t, x, xt, u, ux, uxx, ut, uxt) if index > 0
c
c           for k = 1, ..., npde
c
c DEFIVS:U  subroutine DEFIVS(npde, x, u, ux) 
c           integer          npde
c           double precision x, u(npde), ux(npde)
c 
c           this subroutine is used to define the initial values for 
c           the physical solution and its first spatial derivative, 
c           where u and ux must be defined as functions of x (input 
c           variable)
c
c DEFMSH:U  subroutine DEFMSH(npts, xmesh)
c           integer          npts
c           double precision xmesh(npts)
c
c           this subroutine is used to define the initial mesh, where
c           the entire initial mesh xmesh(1) < xmesh(2) < ...
c           < xmesh(npts) must be defined if job = 1.
c
c DEFBCM:U  subroutine DEFBCM(npde, index, t, x, xt, u, ux, uxx, ut,
c                             uxt, res)
c           integer          npde, index
c           double precision t, x, xt, u(npde), ux(npde), uxx(npde),
c                            ut(npde), uxt(npde), res
c 
c           this subroutine is used to define the boundary conditions
c           for the mesh (the coordinate transformation). the residual
c           res must be defined in terms of input variables npde, index,
c           t, x, xt, u, ux, uxx, ut and uxt by 
c
c           res := b^L_M(t, x, xt, u, ux, uxx, ut, uxt)   if index < 0
c           res := b^R_M(t, x, xt, u, ux, uxx, ut, uxt)   if index > 0
c
c           NOTE: for fixed endpoints, simply use res := x - x^L if 
c           index < 0 and res := x - x^R if index > 0  
c
c DEFMNT:U  subroutine DEFMNT(npde, t, x, u, ux, uxx, ut, uxt, fmntr)
c           integer          npde 
c           double precision t, x, u(npde), ux(npde), uxx(npde),
c                            ut(npde), uxt(npde), fmntr
c
c           this subroutine is used to define the monitor function. the
c           value fmntr of the monitor function must be defined in 
c           terms of input variables npde, t, x, u, ux, uxx, ut and uxt
c
c           NOTE: for the arclength monitor function,
c           fmntr := dsqrt( 1 + ux(1)^2 + ... + ux(npde)^2 )
c
c DEFOUT:U  subroutine DEFOUT(npde, npts, t, xmesh, xmesht, u, ux, ut,
c                             tstep, touta, ntouta, istop, index, nts)
c           integer          npde, npts, i, ntouta, istop, index, nts
c           double precision xmesh(npts), xmesht(npts), u(npts,npde),
c                            ux(npts,npde), ut(npts,npde), tstep,
c                            touta(ntouta)
c
c           this subroutine is used to output the mesh points (xmesh) 
c           and mesh speeds (xmesht) and the solution values (u, ut, ux) 
c           at time t, where u(i,k) := u_k(x_i(t),t) (k=1,...,mpdes;
c           i=1,...,npts), tstep is the current time stepsize and touta 
c           is the array of output times prescribed by the user. the user 
c           can stop the time integration when some condition is met by 
c           setting istop < 0.  c           index is a variable such that 
c
c           the code generates a mesh                   when index < 0
c                               (for job = 2 or 3) 
c           the code is solves the physical pde         when index > 0 
c                               (for job = 1 or 2) 
c
c           all variables are input variables 
c
c DEFFNC:I  subroutine DEFFNC(x, k, npde, npts, xmesh, u, ux, uval,
c                             uxval)
c           integer          k, npde, npts
c           double precision x, xmesh(npts), u(npts,npde),
c                            ux(npts,npde), uval, uxval
c
c           this subroutine provides the value of the k_th component 
c           of the solution and its first spatial derivative at x, 
c           using a cubic Hermite spline which interpolates the
c           given solution values (u and ux) at mesh points (xmesh)
c           (see DEFOUT). it is provided for convenient user 
c           evaluation of the solution and its first spatial 
c           derivative at non-mesh points 
c
c MOVCL1:I  this is the key subroutine (MOVCOL is simply a user
c           interface subroutine)
c
c RESODE:I  defines the ode system for DASSL
c
c JACODE:I  defines a dummy routine for DASSL
c
c RESPDE:I  computes the residuals for the discretizations of the 
c           physical pdes and their bcs
c
c RESPD1:I  computes the residuals for the discretizations of the 
c           (artificial) pde and bcs during the mesh generation
c
c EVLMNT:I  computes the monitor function 
c
c SMTMNT:I  smooths the monitor function (only used for mmpde = 2 
c           and 3)
c
c RESMEQ:I  computes the residuals of the discrete mesh equations and 
c           their bcs 
c
c DRVTVS:I  computes values of u, u_x, u_xx, u_t and u_xt for the 
c           cubic Hermite interpolatory spline solution at point x
c
c FNCSHP:I  defines the shape function and derivative values for 
c           cubic Hermite polynomial interpolation
c
c ZERMCH:I  computes the machine unit roundoff in double precision
c
c SLNOUT:I  driver for the subroutine DEFOUT for solution output 
c
c DDASSL    double precision version of DASSL. this and certain
c           subroutines (from LINPACK) used by DASSL must be included. 
c           for a description see [Petzold]
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                     c
c END OF THE DESCRIPTION OF VARIABLES                                 c
c                                                                     c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c...this subroutine is the code's user interface which checks input 
c...parameters, sets the default values, simplifies the declaration 
c...of workspace arrays and calls the main subroutine movcl1
c
c
      implicit real*8 (a-h, o-z)
      dimension iwork(liw)
      dimension par(20), rwork(lrw), touta(ntouta)
c
      data zero/0.d0/, one/1.d0/, tau1/1.d-4/
c
      zrmch = zermch(dum)
      iflag = 0
c
c...check basic parameters job, inform, mmpde, tau, gamma, ip and
c...stpmax. if their values are not within the permitted ranges,
c...default values are set and iflag = 1. 
c
      job = par(1) + zrmch
      if (par(1).le.-zrmch) job = -1
      if (job.le.0 .or. job.ge.4) then
         iflag = 1
         job = 2
      endif
c
      inform = par(2) + zrmch
      if (par(2).le.-zrmch) inform = -1
      if (inform.le.0) then
         iflag = 1
         inform = 6
      endif
c
      mmpde = par(3) + zrmch
      if (par(3).le.-zrmch) mmpde = -1
      if (mmpde.le.0 .or. mmpde.ge.5) then
         iflag  = 1
         mmpde  = 4
      endif
c
      tau = par(4) + zrmch
      if (tau.le.zero) then
         iflag = 1
         tau = tau1
      endif
c
      gamma = par(5) + zrmch
      if (gamma.le.zero) then
         iflag = 1
         gamma = one
      endif
c
      ip = par(6) + zrmch
      if (par(6).le.-zrmch) ip = -1
      if (ip.lt.0) then
         iflag = 1
         ip = 2
      endif
c
      stpmax = par(7)
      if (stpmax.le.zrmch) then
         iflag = 1
         stpmax = one/zrmch
      endif
c
c...check basic parameters npde, npts, atol, rtol, ntouta,
c...xl and xr
c
      if (npde.le.0) iflag = -1
      if (npts.le.0) iflag = -1
      if (dabs(atol)+dabs(rtol).le.zrmch) iflag = -1
      if (ntouta.le.1) iflag = -1
      if (job.eq.2.or.job.eq.3) then
         xl = par(8)
         xr = par(9)
         if (xl.gt.xr) iflag = -1
      endif
c
c...if iflag = - 1, at least one of the parameters npde, npts, 
c...atol, rtol and ntouta has been incorrectly input, so the 
c...computation is halted
c
      if (iflag.eq.-1) then
         write(inform, 99910) iflag
         write (inform, 99990) npde, npts, job, inform, mmpde, tau,
     +   gamma, ip, atol, rtol, stpmax, ntouta
         return
      endif
c
c...check the dimensions of the workspace arrays
c
      neq = (2*npde+1)*npts
      if (mmpde.eq.2 .or. mmpde.eq.3) then
         lrw0 = 40+12*neq+3*(ip+2)*(2*npde+1)*neq
      else
         lrw0 = 40+12*neq+3*(0+2)*(2*npde+1)*neq
      endif
      lrw1 = 20+6*npde+(3*npde+2)*npts
      if (lrw.lt.2*npde+lrw0+lrw1) iflag = -2
      if (liw.lt.20+npts*(2*npde+1)) iflag = -2
c
c...when iflag = - 2, at least one of the dimensions of the
c...workspace arrays rwork and iwork is incorrectly input,
c...so the computation is stopped
c
      if (iflag.eq.-2) then
         write(inform, 99920) iflag
         write (inform, 99990) npde, npts, job, inform, mmpde, tau,
     +   gamma, ip, atol, rtol, stpmax, ntouta
         return
      endif
c
c...call movcl1
c
      call movcl1(npde, npts, job, inform, iflag, xl, xr,
     1 mmpde, tau, gamma, ip, atol, rtol, stpmax, touta, ntouta,
     2 rwork(1), rwork(neq+1), rwork(2*neq+1), lrw0, iwork, liw,
     3 rwork(2*neq+lrw0+1), lrw1)
c
      return
c
99910 format(
     +/' *****************************************************'
     +/' iflag = ',i3
     +/' halt......at least one of the parameters npde, npts'
     +/' atol, rtol, ntouta, xl, and xr is incorrectly input.'
     +/' they must satisfy :'
     +/' npde >= 1                    npts >= 1'
     +/' |atol| + |rtol| > 0          ntouta >= 2'
     +/' xl <= xr (when job = 2 and 3)'
     +/' *****************************************************'/)
99920 format(
     +/' *****************************************************'
     +/' iflag = ',i3
     +/' halt......at least one of the workspace arrays rwork'
     +/' and iwork has insufficient length. check the values'
     +/' of lrw and liw. also check the values for npde, npts'
     +/' and ip := par(6)'
     +/' *****************************************************'/)
99990 format(
     +/' *****************************************************'
     +/'  npde =',i12,'  npts =',i12,
     +'   job =',i12,' inform=', i10
     +/' mmpde =',i12,'   tau =',e12.5,
     +' gamma =',e12.5,'    ip =',i10
     +/'  atol =',e12.5,'  rtol =', e12.5,
     +' stpmax=',e12.4,' ntouta=',i10
     +/' *****************************************************'/)
      end
c
c***********************************************************************
c***********************************************************************
c
      subroutine movcl1(npde, npts, job, inform, iflag, xl, xr,
     +mmpde, tau, gamma, ip, atol, rtol, stpmax, touta, ntouta,
     +y, ydot, rwk, lrw, iwk, liw, rwk1, lrw1)
c
c...this is the main subroutine. 
c...the relationships between y, ydot, and (x, u) are
c...
c...       y(m*i)               := x_i(t)
c...       y(m*(i-1)+2*(k-1)+1) := u_k(x_i(t),t)
c...       y(m*(i-1)+2*(k-1)+2) := (u_k)_x(x_i(t),t)
c...
c...                               k = 1, ..., npde
c...                               i = 1, ..., npts
c...
c...       ydot := y_t
c...
c...and m = 2*npde+1
c...
c...job = 1: solve the physical pdes starting with the given mesh
c...job = 2: first generate a mesh corresponding to the given
c...         initial value of the solution and then solve
c...         the physical pdes starting with this mesh
c...job = 3: only generate a mesh corresponding to the given
c...         initial value of the solution
c
      implicit real*8 (a-h, o-z)
      dimension iwk(liw), info(15), iwk1(20)
      dimension touta(ntouta), y(npts*(2*npde+1)), rwk(lrw)
      dimension rwk1(lrw1), ydot(npts*(2*npde+1))
      logical phypde
      external  resode, jacode
c
      data zero/0.d0/, half/0.5d0/, one/1.d0/
c
c...define some basic parameters
c
      m = 2*npde+1
      neq = m*npts
      m11 = 20
      m12 = 20 + 1*npde
      m13 = 20 + 2*npde
      m14 = 20 + 3*npde
      m15 = 20 + 4*npde
      m16 = 20 + 5*npde
      m21 = 20 + 6*npde
      m22 = 20 + 6*npde + 1*npts
      m31 = 20 + 6*npde + 2*npts
      m32 = 20 + 6*npde + 2*npts + 1*npde*npts
      m33 = 20 + 6*npde + 2*npts + 2*npde*npts
      zrmch = zermch(dum)
c     phypde = 'true' for solving the physical pdes
c     phypde = 'false' for generating a mesh
      if (job.eq.1) phypde = .true.
      if (job.eq.2) phypde = .false.
      if (job.eq.3) phypde = .false.
c
c...compute the initial mesh
c
c     for solving the physical pde
      if (phypde) then
         call defmsh(npts, rwk1(m21+1))
         do 10 i = 1, npts
            y(m*i) = rwk1(m21+i)
   10    continue
      endif
c
c     the mesh generation case: starting from a uniform mesh
      if (.not.phypde) then
         do 15 i = 1, npts
            y(m*i) = xl + (xr-xl)*dfloat(i-1)/dfloat(npts-1)
   15    continue
      endif
c
c...print the headers
c
      write (inform, 99920)
   20 if (phypde) then
         write(inform, 99930)
      else
         write(inform, 99940)
      endif
c
c...values for t, y and ydot which are consistent with the ODE
c...system (arising in the method of lines procedure) are 
c...needed initially to ensure convergence of ddassl at the first 
c...step of the numerical integration. these initial values are 
c...first obtained 
c
c...the initial time derivatives ydot are computed by taking two 
c...time integration steps with the first-order backward 
c...differentiation formula and suppressing ddassl's time 
c...truncation error control during these steps
c
c...compute the initial values of the physical solution
c
      if (phypde) then
c        the case of solving the physical pdes
         do 25 i = 1, npts
         call defivs(npde, y(m*i), rwk1(m11+1), rwk1(m12+1))
         do 25 k = 1, npde
         y(m*(i-1)+2*(k-1)+1) = rwk1(m11+k)
         y(m*(i-1)+2*(k-1)+2) = rwk1(m12+k)
   25    continue
      else
c        the mesh generation case
         do 30 i = 1, npts
         do 30 k = 1, npde
         y(m*(i-1)+2*(k-1)+1) = zero 
         y(m*(i-1)+2*(k-1)+2) = zero 
   30    continue
      endif
c
c...define the input parameters for subroutine resode
c
      iwk1(1) = npde
      iwk1(2) = npts
      iwk1(3) = mmpde
      iwk1(4) = ip
      if (phypde) then
         iwk1(5) = + 1
      else
         iwk1(5) = - 1
      endif
      rwk1(1) = gamma
      rwk1(2) = tau
c
c...define the input parameters for ddassl
c
      do 35 i = 1, 15
      info(i) = 0
   35 continue
      info(3) = 1
      info(6) = 1
c        define the lower and upper bandwidths for the jacobian
         if (mmpde.eq.2 .or. mmpde.eq.3) then
            iwk(1) = min0(m*(ip+2), neq-1)
         else
            iwk(1) = min0(m*(0+2), neq-1)
         endif
         iwk(2) = iwk(1)
      info(7) = 1
c        define the maximal time stepsize
         rwk(2) = dmin1(dsqrt(zrmch), atol, rtol)
      info(8) = 1
c        define the initial time stepsize
         rwk(3) = rwk(2)*half
      info(9) = 1
c        define the order of the bdf method (1st order)
         iwk(3) = 1
      info(11) = 1
c        define the initial values for ydot
         do 40 i = 1, neq
            ydot(i) = zero 
   40    continue
c
c...take two time integration steps using ddassl
c
      if (phypde) then
         t = touta(1)
         tout = touta(1)+one
      else
         t = zero 
         tout = one 
      endif
      atol1 = one 
      rtol1 = one 
      do 50 itout = 1, 2
      call ddassl(
     +resode, neq, t, y, ydot, tout, info, rtol1, atol1,
     +idid, rwk, lrw, iwk, liw, rwk1, iwk1, jacode)
   50 continue
c
c...now begin the actual numerical integration
c
c...the physical and mesh pdes are discretized in space with
c...a moving collocation method. the resulting ode system is
c...integrated in time with ddassl
c
c...reset the initial data for the mesh and solution
c
c     for the mesh
c
c     for solving the physical pde
      if (job.eq.1) then
         call defmsh(npts, rwk1(m21+1))
         do 54 i = 1, npts
            y(m*i) = rwk1(m21+i)
   54    continue
      endif
c
c     the mesh generation case: starting from a uniform mesh
      if (.not.phypde) then
         do 55 i = 1, npts
            y(m*i) = xl + (xr-xl)*dfloat(i-1)/dfloat(npts-1)
   55    continue
      endif
c
c     for the physical solution
      if (phypde) then
c        the case of solving the physical pdes
         do 56 i = 1, npts
         call defivs(npde, y(m*i), rwk1(m11+1), rwk1(m12+1))
         do 56 k = 1, npde
         y(m*(i-1)+2*(k-1)+1) = rwk1(m11+k)
         y(m*(i-1)+2*(k-1)+2) = rwk1(m12+k)
   56    continue
      else
c        the mesh generation case
         do 60 i = 1, npts
         do 60 k = 1, npde
         y(m*(i-1)+2*(k-1)+1) = zero 
         y(m*(i-1)+2*(k-1)+2) = zero 
   60    continue
      end if
c
c...define the input parameters for subroutine resode
c
      iwk1(1) = npde
      iwk1(2) = npts
      iwk1(3) = mmpde
      iwk1(4) = ip
      if (phypde) then
         iwk1(5) = + 1
      else
         iwk1(5) = - 1
      endif
      rwk1(1) = gamma
      rwk1(2) = tau
c
c...define the input parameters for ddassl
c
      do 65 i = 1, 15
      info(i) = 0
   65 continue
      info(3) = 1
      info(6) = 1
c        define the lower and upper bandwidths for the jacobian
         if (mmpde.eq.2 .or. mmpde.eq.3) then
            iwk(1) = min0(m*(ip+2), neq-1)
         else
            iwk(1) = min0(m*(0+2), neq-1)
         endif
         iwk(2) = iwk(1)
      info(7) = 1
c        define the maximal time stepsize
         rwk(2) = stpmax
      info(8) = 1
c        define the initial time stepsize
         temp = dabs(touta(ntouta))
         if (.not.phypde) temp = one 
         rwk(3) = dmin1(stpmax, rtol, atol, temp, one)
      info(11) = 1
c
c...print the initial solution
c
      if (phypde) then
         t = touta(1)
         index = + 1
      else
         t = zero
         index = - 1
      endif
      tstep = zero
      istop = 0
      nts = 0
      call slnout(npde, npts, t, y, ydot, touta, ntouta,
     +            tstep, istop, index, nts, rwk1(m11+1), rwk1(m12+1),
     +            rwk1(m13+1), rwk1(m14+1), rwk1(m15+1),
     +            rwk1(m21+1), rwk1(m22+1), rwk1(m31+1), rwk1(m32+1),
     +            rwk1(m33+1))
c     if istop < 0, stop the computation
      if (istop.lt.0) then
         iflag = -3
         write(inform, 99950) iflag
         write (inform, 99990) npde, npts, job, inform, mmpde, tau,
     +   gamma, ip, atol, rtol, stpmax, ntouta
         return
      endif
c
c...perform the time integration with ddassl
c
      if (phypde) then
         itoutm = ntouta
      else
         itoutm = 2
      endif
      do 130 itout = 2, itoutm
      if (phypde) then
         tout = touta(itout)
      else
         tout = one
      endif
  110 call ddassl(
     +resode, neq, t, y, ydot, tout, info, rtol, atol,
     +idid, rwk, lrw, iwk, liw, rwk1, iwk1, jacode)
c
c...print the solution at t
c
      if (phypde) then
         index = + 1
      else
         index = - 1
      endif
      tstep = rwk(7)  
      istop = 0
      nts = iwk(11)
      call slnout(npde, npts, t, y, ydot, touta, ntouta,
     +            tstep, istop, index, nts, rwk1(m11+1), rwk1(m12+1),
     +            rwk1(m13+1), rwk1(m14+1), rwk1(m15+1),
     +            rwk1(m21+1), rwk1(m22+1), rwk1(m31+1), rwk1(m32+1),
     +            rwk1(m33+1))
c     if istop < 0, stop the computation
      if (istop.lt.0) then
         iflag = -3
         write(inform, 99950) iflag
         write (inform, 99990) npde, npts, job, inform, mmpde, tau,
     +   gamma, ip, atol, rtol, stpmax, ntouta
         return
      endif
c
c...check the performance of the time integrator
c
      if(idid.eq.1) go to 110
c        a step was successfully taken in the intermediate-output
c        mode. ddassl has not yet reached tout. continue the time
c        integration
      if (idid.eq.-1) then
c        a large amount of work has been expended. reset the
c        limitation and continue the time integration
         info(1) = 1
         goto 110
      endif
      if (idid.ge.-1 .and. idid.le.3) then
c        the time integration is successful to tout. output the
c        performance information
         write (inform, 99970) t, idid, iwk(14), iwk(15),
     +   iwk(11), iwk(13), iwk(12), iwk(8), rwk(7), rwk(3)
         info(1) = 1
      else
c        the time integration fails and the computation is stopped
         iflag = -4
         write (inform, 99980) iflag, idid
         write (inform, 99990) npde, npts, job, inform, mmpde, tau,
     +   gamma, ip, atol, rtol, stpmax, ntouta
         return
      endif
  130 continue
c
c...if job = 2, go back to solve the physical pdes
c
      if (job.eq.2 .and. (.not.phypde)) then
         phypde = .true.
         go to 20
      endif
c
c...the computation is completed. output summary information and
c...return
c
      write (inform, 99990) npde, npts, job, inform, mmpde, tau,
     +gamma, ip, atol, rtol, stpmax, ntouta
      return
99920 format(
     +/' *****************************************************'
     +/' ***                                               ***'
     +/' ***                  MOVCOL                       ***'
     +/' ***                                               ***'
     +/' *****************************************************'
     +/' NTS  : time steps taken'
     +/' JAC  : Jacobian evaluations'
     +/' ETF  : error test failures'
     +/' CFN  : convergence test failures in Newton iteration'
     +/' IDID : indicator of the performance of DASSL'
     +/' NRES : number of the ODE residual evaluations'
     +/' ORD  : order of the method for the time integration'
     +/' STEPLAST : time step used in last step'
     +/' STEPNOW :  time step used in current step'
     +/' *****************************************************'/)
99930 format(
     +/' *****************************************************'
     +/' ***                                               ***'
     +/' ***           solve the physical PDEs             ***'
     +/' ***                                               ***'
     +/' *****************************************************'/)
99940 format(
     +/' *****************************************************'
     +/' ***                                               ***'
     +/' ***       generate an equidistributed mesh        ***'
     +/' ***                                               ***'
     +/' *****************************************************'/)
99950 format(
     +/' *****************************************************'
     +/' iflag = ',i3
     +/' halt......istop is set to be < 0 by the user in DEFOUT'
     +/' *****************************************************'/)
99970 format(
     +/' ---------------------- information for time integration',
     +' ---------------------'/
     +'  time =',e12.5,'  IDID =',i12,'   ETF =',i12,'   CFN =',i10/
     +'   NTS =',i12,'   JAC =',i12,'  NRES =',i12,'   ORD =',i10/
     +' STEPLAST=',e12.5,' STEPNOW=',e12.5,
     +/' -----------------------------------------------------',
     +'------------------------')
99980 format(
     +/' *****************************************************'
     +/' iflag = ',i3,' IDID =', i3
     +/' halt......IDID <= -2 in DASSL (see documentation in'
     +/' DASSL)'
     +/' *****************************************************'/)
99990 format(
     +/' *****************************************************',
     +'************************'
     +/'  npde =',i12,'  npts =',i12,
     +'   job =',i12,' inform=', i10
     +/' mmpde =',i12,'   tau =',e12.5,
     +' gamma =',e12.5,'    ip =',i10
     +/'  atol =',e12.5,'  rtol =', e12.5,
     +' stpmax=',e12.4,' ntouta=',i10
     +/' *****************************************************',
     +'************************'/)
      end
c
c***********************************************************************
c***********************************************************************
c
      subroutine jacode(t, y, ydot, pd, cj, rwk, iwk)
c
c...this dummy routine is required by ddassl
c
      implicit real*8 (a-h, o-z)
      dimension y(*), ydot(*), pd(10,*), rwk(*)
      dimension iwk(*)
      return
      end
c
c***********************************************************************
c***********************************************************************
c
      subroutine resode(t, y, ydot, res, ires, rwk, iwk)
c
c...this subroutine used by ddassl computes the residuals for 
c...the ode system resulting from the discretization of the 
c...physical and mesh pdes. 
c...the relationships between y, ydot, res c...and (x, u) are
c...
c...       y(m*i)               := x_i(t)
c...       y(m*(i-1)+2*(k-1)+1) := u_k(x_i(t),t)
c...       y(m*(i-1)+2*(k-1)+2) := (u_k)_x(x_i(t),t)
c...
c...                               k = 1, ..., npde
c...                               i = 1, ..., npts
c...
c...       ydot := y_t
c...       res  := F(t, x, u, ux, ut, uxt)-[ G(t, x, u, ux, ut, uxt) ]_x
c...
c...and m = 2*npde+1
c
      implicit real*8 (a-h, o-z)
      dimension iwk(*)
      dimension y(*), ydot(*), res(*), rwk(*)
      logical phypde
c
c...define some basic parameters
c
      npde = iwk(1)
      npts = iwk(2)
      mmpde = iwk(3)
      ip = iwk(4)
      if (iwk(5).gt.0) then
c        the case for solving the physical pdes
         phypde = .true.
      else
c        the mesh generation case
         phypde = .false.
      endif
      gamma = rwk(1)
      tau = rwk(2)

      m = 2*npde+1
      m11 = 20
      m12 = 20+1*npde
      m13 = 20+2*npde
      m14 = 20+3*npde
      m15 = 20+4*npde
      m16 = 20+5*npde
      m21 = 20+6*npde
      m22 = 20+6*npde+npts
c
c...compute residuals of the physical pdes and their bcs
c
      if (phypde) then 
c        the case of solving the physical pdes
         call respde(npde, npts, t, y, ydot, res, rwk(m11+1),
     +               rwk(m12+1), rwk(m13+1), rwk(m14+1), rwk(m15+1),
     +               rwk(m16+1))
      else
c        the mesh generation case
         call respd1(npde, npts, t, y, ydot, res, rwk(m11+1),
     +               rwk(m12+1), rwk(m13+1), rwk(m14+1), rwk(m15+1))
      endif
c
c...compute and smooth the monitor function
c
c     rwk(m21+i) := fmntr(i)
      if (mmpde.ge.2)
     +call evlmnt(npde, npts, t, y, ydot, rwk(m21+1), rwk(m11+1),
     +            rwk(m12+1), rwk(m13+1), rwk(m14+1), rwk(m15+1))
      if (mmpde.eq.2 .or. mmpde.eq.3) 
     +call smtmnt(npts, ip, rwk(m21+1), rwk(m22+1))
c
c...compute residuals for the discrete mesh equations and their bcs.
c
      call resmeq(npde, npts, mmpde, tau, gamma, y, ydot, res,
     +            rwk(m21+1), rwk(m11+1), rwk(m12+1), rwk(m13+1),
     +            rwk(m14+1), rwk(m15+1))
      if (.not.phypde) then 
c        the mesh generation case (for which boundaries are fixed)
         res(m) = ydot(m)
         res(m*npts) = ydot(m*npts)
      endif
c
      return
      end
c
c***********************************************************************
c***********************************************************************
c
      subroutine respde(npde, npts, t, y, ydot, res,
     +                  u, ux, uxx, ut, uxt, rw)
c
c...this subroutine computes the residuals for the discretizations
c...of the physical pdes and their bcs. 
c...the relationships between y, ydot, res, and (x, u) are
c...
c...       y(m*i)               := x_i(t)
c...       y(m*(i-1)+2*(k-1)+1) := u_k(x_i(t),t) 
c...       y(m*(i-1)+2*(k-1)+2) := (u_k)_x(x_i(t),t)
c...
c...                               k = 1, ..., npde
c...                               i = 1, ..., npts
c...
c...       ydot := y_t 
c...       res  := F(t, x, u, ux, ut, uxt)-[ G(t, x, u, ux, ut, uxt) ]_x
c...
c...and m = 2*npde+1
c
      implicit real*8 (a-h, o-z)
      dimension y((2*npde+1)*npts), ydot((2*npde+1)*npts)
      dimension res((2*npde+1)*npts), u(npde), ux(npde), uxx(npde)
      dimension ut(npde), uxt(npde), rw(npde), w(3,2)
c
      data s1/0.211324865347452d0/, s2/0.788675134652548d0/
      data half/0.5d0/, one/1.d0/, two/2.d0/, four/4.d0/
c
      m = 2*npde+1
c
c...define the weights for the cell-average collocation approximation
c...to the right-hand-side term
c
      w(1,1) = -one-two*(s2-s1)
      w(2,1) = four*(s2-s1) 
      w(3,1) = one-two*(s2-s1)
      w(1,2) = -one+two*(s2-s1)
      w(2,2) = -four*(s2-s1) 
      w(3,2) = one+two*(s2-s1)
c
c...compute the residuals for the discretization of the physical pdes 
c...at the interior (Gauss) points
c
      do 50 i = 1, npts-1
      h = y(m*(i+1))-y(m*i)
c
c...compute the left-hand-side term F(...) using collocation
c
      do 20 j = 1, 2
      if (j.eq.1) x = y(m*i)+s1*h
      if (j.eq.2) x = y(m*i)+s2*h
      call drvtvs(npde, npts, i, x, y, ydot, u, ux, uxx, ut, uxt)
      call defpde(npde, -1, t, x, u, ux, ut, uxt, rw)
      do 10 k = 1, npde
      if(j.eq.1) res(m*(i-1)+2*(k-1) +2) = rw(k)
      if(j.eq.2) res(m*(i  )+2*(k-1) +1) = rw(k)
   10 continue
   20 continue
c
c...compute the right-hand-side term [ G(...) ]_x using cell-average
c...collocation
c
      do 40 j = 1, 3
      if (j.eq.1) x = y(m*i) 
      if (j.eq.2) x = y(m*i)+half*h
      if (j.eq.3) x = y(m*(i+1))
      call drvtvs(npde, npts, i, x, y, ydot, u, ux, uxx, ut, uxt)
      call defpde(npde, +1, t, x, u, ux, ut, uxt, rw)
      do 30 k = 1, npde
      res(m*(i-1)+2*(k-1) +2) = res(m*(i-1)+2*(k-1) +2)-rw(k)*w(j,1)/h
      res(m*(i  )+2*(k-1) +1) = res(m*(i  )+2*(k-1) +1)-rw(k)*w(j,2)/h
   30 continue
   40 continue
   50 continue
c
c...compute the bcs at x = x^L (left end) using collocation
c
      x = y(m)
      xt = ydot(m)
      call drvtvs(npde, npts, 1, x, y, ydot, u, ux, uxx, ut, uxt)
      call defbcp(npde, -1, t, x, xt, u, ux, uxx, ut, uxt, rw)
      do 60 k = 1, npde
      res(2*(k-1)+1) = rw(k)
   60 continue
c
c...compute the bcs at x = x^R (right end) using collocation
c
      x = y(m*npts)
      xt = ydot(m*npts)
      call drvtvs(npde, npts, npts-1, x, y, ydot, u, ux, uxx, ut,
     +uxt)
      call defbcp(npde, +1, t, x, xt, u, ux, uxx, ut, uxt, rw)
      do 70 k = 1, npde
      res(m*(npts-1)+2*(k-1)+2) = rw(k)
   70 continue
c
      return
      end
c
c***********************************************************************
c***********************************************************************
c
      subroutine respd1(npde, npts, t, y, ydot, res,
     +                  u, ux, uxx, ut, uxt)
c
c...for the mesh generation case, the (artificial) physical pde is
c...defined as
c...          u_t = U(x), 0 < t <= 1
c...
c...where U(x) is the given initial solution defined in subroutine
c...defivs
c... 
c...this subroutine computes the residuals for the discretization 
c...to this pde and for the bcs.
c...the relationships between y, ydot, res, and (x, u) are
c...
c...       y(m*i)               := x_i(t)
c...       y(m*(i-1)+2*(k-1)+1) := u_k(x_i(t),t) 
c...       y(m*(i-1)+2*(k-1)+2) := (u_k)_x(x_i(t),t)
c...
c...                               k = 1, ..., npde
c...                               i = 1, ..., npts
c...
c...       ydot := y_t 
c...       res  := u_t - U(x) 
c...
c...and m = 2*npde+1
c
      implicit real*8 (a-h, o-z)
      dimension y((2*npde+1)*npts), ydot((2*npde+1)*npts)
      dimension res((2*npde+1)*npts), u(npde), ux(npde), uxx(npde)
      dimension ut(npde), uxt(npde)
c
      data s1/0.211324865347452d0/, s2/0.788675134652548d0/
      data zero/0.d0/, one/1.d0/
c
      m = 2*npde+1
c
c...compute the residuals for the discretizations to the pdes at 
c...the interior (Gauss) points
c
      do 30 i = 1, npts-1
      h = y(m*(i+1))-y(m*i)
      ht = ydot(m*(i+1))-ydot(m*i)
      do 20 j = 1, 2
      if (j.eq.1) then
         x = y(m*i)+s1*h
         xt = ydot(m*i)+s1*ht
      endif
      if (j.eq.2) then
         x = y(m*i)+s2*h
         xt = ydot(m*i)+s2*ht
      endif
      call drvtvs(npde, npts, i, x, y, ydot, u, ux, uxx, ut, uxt)
      call defivs(npde, x, u, uxx)
      do 10 k = 1, npde
c
c...1st order upwind treatment
c
      ut(k) = ut(k) + ux(k)*xt
      temp0 = y(m*(i-1)+2*(k-1)+1)*fncshp(1,0,zero)
     +       +y(m*(i-1)+2*(k-1)+2)*fncshp(2,0,zero)*h
     +       +y(m*(i  )+2*(k-1)+1)*fncshp(3,0,zero)
     +       +y(m*(i  )+2*(k-1)+2)*fncshp(4,0,zero)*h
      temp1 = y(m*(i-1)+2*(k-1)+1)*fncshp(1,0,s1)
     +       +y(m*(i-1)+2*(k-1)+2)*fncshp(2,0,s1)*h
     +       +y(m*(i  )+2*(k-1)+1)*fncshp(3,0,s1)
     +       +y(m*(i  )+2*(k-1)+2)*fncshp(4,0,s1)*h
      temp2 = y(m*(i-1)+2*(k-1)+1)*fncshp(1,0,s2)
     +       +y(m*(i-1)+2*(k-1)+2)*fncshp(2,0,s2)*h
     +       +y(m*(i  )+2*(k-1)+1)*fncshp(3,0,s2)
     +       +y(m*(i  )+2*(k-1)+2)*fncshp(4,0,s2)*h
      temp3 = y(m*(i-1)+2*(k-1)+1)*fncshp(1,0,one)
     +       +y(m*(i-1)+2*(k-1)+2)*fncshp(2,0,one)*h
     +       +y(m*(i  )+2*(k-1)+1)*fncshp(3,0,one)
     +       +y(m*(i  )+2*(k-1)+2)*fncshp(4,0,one)*h
      if (j.eq.1) then
         if (xt.ge.zero) temp = (temp2-temp1)/(s2-s1)/h
         if (xt.le.zero) temp = (temp1-temp0)/s1/h
         res(m*(i-1)+2*(k-1) +2) = (ut(k)-temp*xt) -u(k) 
      endif
      if (j.eq.2) then
         if (xt.ge.zero) temp = (temp3-temp2)/(one-s2)/h
         if (xt.le.zero) temp = (temp2-temp1)/(s2-s1)/h
         res(m*(i  )+2*(k-1) +1) = (ut(k)-temp*xt) -u(k) 
      endif
   10 continue
   20 continue
   30 continue
c
c...compute the residual for the collocation solution of the
c...pdes at x = x^L (left end) 
c
      x  = y(m)
      xt = ydot(m)
      call drvtvs(npde, npts, 1, x, y, ydot, u, ux, uxx, ut, uxt)
      call defivs(npde, x, ux, uxx)
c     here ux represents U(x)
      do 40 k = 1, npde
      res(2*(k-1)+1) = u(k)-t*ux(k)
   40 continue
c
c...compute the residual for the collocation solution of the
c...pdes at x = x^R (right end) 
c
c
      x  = y(m*npts)
      xt = ydot(m*npts)
      call drvtvs(npde, npts, npts-1, x, y, ydot, u, ux, uxx, ut, uxt)
      call defivs(npde, x, ux, uxx)
c     here ux represents U(x)
      do 50 k = 1, npde
      res(m*(npts-1)+2*(k-1)+2) = u(k)-t*ux(k)
   50 continue
c
      return
      end
c
c***********************************************************************
c***********************************************************************
c
      subroutine resmeq(npde, npts, mmpde, tau, gamma, y, ydot, res,
     +                  fmntr, u, ux, uxx, ut, uxt)
c
c...this subroutine computes the residuals for the finite difference
c...approximations of the mesh equations c...and their bcs. 
c...the relationships between y, ydot, res, and (x, u) are
c...
c...       y(m*i)               := x_i(t)
c...       y(m*(i-1)+2*(k-1)+1) := u_k(x_i(t),t) 
c...       y(m*(i-1)+2*(k-1)+2) := (u_k)_x(x_i(t),t)
c...
c...                               k = 1, ..., npde
c...                               i = 1, ..., npts
c...
c...       ydot := y_t 
c...
c...and m = 2*npde+1
c
      implicit real*8 (a-h, o-z)
      dimension y((2*npde+1)*npts), ydot((2*npde+1)*npts)
      dimension res((2*npde+1)*npts), fmntr(npts), u(npde)
      dimension ux(npde), uxx(npde), ut(npde), uxt(npde) 
c
      data one/1.d0/, two/2.d0/
c
      m = 2*npde+1
c
c...mmpde = 1: fixed mesh
c
      if (mmpde.eq.1) then
         do 10 i = 1, npts
         res(m*i) = ydot(m*i)
   10    continue
         return
      endif
c
c...compute bcs for the coordinate transformation (the mesh)
c...for mpde = 2, 3, and 4
c
      x = y(m)
      xt = ydot(m)
      call drvtvs(npde, npts, 1, x, y, ydot, u, ux, uxx, ut, uxt)
      call defbcm(npde, -1, t, x, xt, u, ux, uxx, ut, uxt, temp)
      res(m) = temp 
      x = y(m*npts)
      xt = ydot(m*npts)
      call drvtvs(npde, npts, npts-1, x, y, ydot, u, ux, uxx, ut, uxt)
      call defbcm(npde, +1, t, x, xt, u, ux, uxx, ut, uxt, temp)
      res(m*npts) = temp
c
c...mmpde = 2: MMPDE4 in [HRR1, HRR2] at the interior mesh points 
c
      if (mmpde.eq.2) then
         do 20 i = 2, npts-1
         temp1 = fmntr(i-1)
         temp2 = fmntr(i)
         res(m*i) =
     +    temp2*(y(m*(i+1))-y(m*i))-temp1*(y(m*i)-y(m*(i-1)))
     +   +( temp2*(ydot(m*(i+1))-ydot(m*i))
     +     -temp1*(ydot(m*i)-ydot(m*(i-1))))*tau
   20    continue
         return
      endif
c
c...mmpde = 3: MMPDE6 in [HRR1, HRR2] at the interior mesh points 
c
      if (mmpde.eq.3) then
         do 30 i = 2, npts-1
         temp1 = fmntr(i-1)
         temp2 = fmntr(i)
         res(m*i) =
     +    temp2*(y(m*(i+1))-y(m*i))-temp1*(y(m*i)-y(m*(i-1)))
     +   +( (ydot(m*(i+1))-ydot(m*i))
     +     -(ydot(m*i)-ydot(m*(i-1))))*tau
   30    continue
         return
      endif
c
c...mmpde = 4: MMPDE(24) in [HR2], which is a smoothed version 
c...of MMPDE4 in terms of the mesh concentration function at the 
c...interior mesh points 
c
      xlambd = gamma*(gamma+one)
      if (mmpde.eq.4) then
         do 40 i = 2, npts-1
         temp1  = fmntr(i-1)
         temp2  = fmntr(i)
         ym1 = -(ydot(m*i)-ydot(m*(i-1)))/(y(m*i)-y(m*(i-1)))**2
     +         *tau+one/(y(m*i)-y(m*(i-1)))
         yp1 = -(ydot(m*(i+1))-ydot(m*i))/(y(m*(i+1))-y(m*i))**2
     +         *tau+one/(y(m*(i+1))-y(m*i))
         if (i-2.ge.1) then
            ym3 = -(ydot(m*(i-1))-ydot(m*(i-2)))
     +             /(y(m*(i-1))-y(m*(i-2)))**2*tau
     +            +one/(y(m*(i-1))-y(m*(i-2)))
         else
            ym3  = ym1
         endif
         if (i+2.le.npts) then
            yp3 = -(ydot(m*(i+2))-ydot(m*(i+1)))
     +             /(y(m*(i+2))-y(m*(i+1)))**2*tau
     +            +one/(y(m*(i+2))-y(m*(i+1)))
         else
            yp3 = yp1
         endif
         res(m*i) = (yp1-xlambd*(yp3-two*yp1+ym1))/temp2
     +                -(ym1-xlambd*(yp1-two*ym1+ym3))/temp1
   40    continue
         return
      endif
c
      return
      end
c
c***********************************************************************
c***********************************************************************
c
      subroutine drvtvs(npde, npts, i, x, y, ydot, u, ux, uxx, ut, uxt)
c
c...this subroutine computes the value (u) and partial derivatives (ux,
c...uxx, ut, uxt) of the cubic Hermite interpolatory spline at any
c...point (x) in the i_th interval [x_i, x_(i+1)], i = 1, ..., npts-1,
c...by using nodal values y and ydot which are defined as
c...
c...       y(m*i)               := x_i(t)
c...       y(m*(i-1)+2*(k-1)+1) := u_k(x_i(t),t) 
c...       y(m*(i-1)+2*(k-1)+2) := (u_k)_x(x_i(t),t)
c...
c...                               k = 1, ..., npde
c...                               i = 1, ..., npts
c...
c...       ydot := y_t 
c...
c...and m = 2*npde+1
c
      implicit real*8 (a-h, o-z)
      dimension y((2*npde+1)*npts), ydot((2*npde+1)*npts)
      dimension u(npde), ux(npde), uxx(npde), ut(npde), uxt(npde)
      dimension temp(4, 0:2)
c
      m = 2*npde+1
      h  = y(m*(i+1))-y(m*i)
      ht = ydot(m*(i+1))-ydot(m*i)
      s = (x-y(m*i))/h
c
      do 5 j = 1, 4
      do 5 k = 0, 2
         temp(j,k) = fncshp(j,k,s)
    5 continue
c
      do 10 k = 1, npde
      u(k) =
     + y(m*(i-1)+2*(k-1)+1)*temp(1,0)
     ++y(m*(i-1)+2*(k-1)+2)*temp(2,0)*h
     ++y(m*(i  )+2*(k-1)+1)*temp(3,0)
     ++y(m*(i  )+2*(k-1)+2)*temp(4,0)*h
      ux(k) = (
     + y(m*(i-1)+2*(k-1)+1)*temp(1,1)
     ++y(m*(i-1)+2*(k-1)+2)*temp(2,1)*h
     ++y(m*(i  )+2*(k-1)+1)*temp(3,1)
     ++y(m*(i  )+2*(k-1)+2)*temp(4,1)*h)/h
      ut(k) =
     + ydot(m*(i-1)+2*(k-1)+1)*temp(1,0)
     ++ydot(m*(i-1)+2*(k-1)+2)*temp(2,0)*h
     ++ydot(m*(i  )+2*(k-1)+1)*temp(3,0)
     ++ydot(m*(i  )+2*(k-1)+2)*temp(4,0)*h
     ++y(m*(i-1)+2*(k-1)+2)*temp(2,0)*ht
     ++y(m*(i  )+2*(k-1)+2)*temp(4,0)*ht
     +-ux(k)*(ydot(m*(i-1)+m)+s*ht)
      uxx(k) = (
     + y(m*(i-1)+2*(k-1)+1)*temp(1,2)
     ++y(m*(i-1)+2*(k-1)+2)*temp(2,2)*h
     ++y(m*(i  )+2*(k-1)+1)*temp(3,2)
     ++y(m*(i  )+2*(k-1)+2)*temp(4,2)*h)/h/h
      uxt(k) =
     +(ydot(m*(i-1)+2*(k-1)+1)*temp(1,1)
     ++ydot(m*(i-1)+2*(k-1)+2)*temp(2,1)*h
     ++ydot(m*(i  )+2*(k-1)+1)*temp(3,1)
     ++ydot(m*(i  )+2*(k-1)+2)*temp(4,1)*h
     ++y(m*(i-1)+2*(k-1)+2)*temp(2,1)*ht
     ++y(m*(i  )+2*(k-1)+2)*temp(4,1)*ht)/h
     +-ux(k)*ht/h
     +-uxx(k)*(ydot(m*(i-1)+m)+s*ht)
   10 continue
c
      return
      end
c
c***********************************************************************
c***********************************************************************
c
      real*8 function fncshp(j, k, s)
c
c...this function subroutine computes, at a point s in [0, 1], the 
c...value for the k_th derivative (k =  0, 1 or 2) of the j_th shape 
c...function (j = 1, 2, 3 or 4) for the cubic Hermite interpolate
c
      implicit real*8 (a-h, o-z)
c
      data one/1.d0/, two/2.d0/, three/3.d0/
      data four/4.d0/, six/6.d0/
c
      if (j.eq.1) then
         if (k.eq.0) fncshp = (one+two*s)*(one-s)*(one-s)
         if (k.eq.1) fncshp = -six*s*(one-s)
         if (k.eq.2) fncshp = -six*(one-two*s)
      endif
      if (j.eq.2) then
         if (k.eq.0) fncshp = s*(one-s)*(one-s)
         if (k.eq.1) fncshp = (one-s)*(one-three*s)
         if (k.eq.2) fncshp = six*s-four 
      endif
      if (j.eq.3) then
         if (k.eq.0) fncshp = (three-two*s)*s*s
         if (k.eq.1) fncshp = six*s*(one-s) 
         if (k.eq.2) fncshp = six*(one-two*s)
      endif
      if (j.eq.4) then
         if (k.eq.0) fncshp = (s-one)*s*s
         if (k.eq.1) fncshp = s*(three*s-two) 
         if (k.eq.2) fncshp = six*s-two
      endif
c
      return
      end
c
c***********************************************************************
c***********************************************************************
c
      subroutine evlmnt(npde, npts, t, y, ydot, fmntr,
     +                  u, ux, uxx, ut, uxt)
c
c...this subroutine computes the monitor function value
c...fmntr(i) := c...M_(i+1/2). 
c...the relationships between y, ydot, and (x, u) are
c...
c...       y(m*i)               := x_i(t)
c...       y(m*(i-1)+2*(k-1)+1) := u_k(x_i(t),t) 
c...       y(m*(i-1)+2*(k-1)+2) := (u_k)_x(x_i(t),t)
c...
c...                               k = 1, ..., npde
c...                               i = 1, ..., npts
c...
c...       ydot := y_t 
c...
c...and m = 2*npde+1
c
      implicit real*8 (a-h, o-z)
      dimension y((2*npde+1)*npts), ydot((2*npde+1)*npts)
      dimension u(npde), ux(npde), uxx(npde), ut(npde), uxt(npde)
      dimension fmntr(npts)
c
      data s1/0.211324865347452d0/, s2/0.788675134652548d0/
      data zero/0.d0/, quart/0.25d0/, half/0.5d0/
c
      m = 2*npde+1
c
      do 20 i = 1, npts-1
      fmntr(i) = zero
      h = y(m*(i+1))-y(m*i)
      do 10 j = 1, 2
      if (j.eq.1) x = y(m*i)+s1*h
      if (j.eq.2) x = y(m*i)+s2*h
      call drvtvs(npde, npts, i, x, y, ydot, u, ux, uxx, ut, uxt)
      call defmnt(npde, t, x, u, ux, uxx, ut, uxt, temp)
      fmntr(i) = fmntr(i) + half*temp
   10 continue
   20 continue
c
c...approximate the monitor function more accurately in the end
c...subintervals (i.e., for i = 1 and npts-1)
c
      do 40 i = 1, npts-1, npts-2
      fmntr(i) = half*fmntr(i)
      do 30 j = i, i+1
      x = y(m*j)
      call drvtvs(npde, npts, i, x, y, ydot, u, ux, uxx, ut, uxt)
      call defmnt(npde, t, x, u, ux, uxx, ut, uxt, temp)
      fmntr(i) = fmntr(i) + quart*temp
   30 continue
   40 continue
c
      return
      end      
c
c***********************************************************************
c***********************************************************************
c
      subroutine smtmnt(npts, ip, fmntr, rw)
c
c...this subroutine smooths the monitor function. it is called 
c...only when mmpde = 2 or mmpde = 3
c
      implicit real*8 (a-h, o-z)
      dimension fmntr(npts), rw(npts)
c
      data gamma/2.d0/, zero/0.d0/, one/1.d0/
c
      do 10 i = 1, npts-1
      rw(i) = fmntr(i)
   10 continue
      do 30 i = 1, npts-1
      temp = zero 
      fmntr(i) = zero 
      do 20 j = max0(1, i-ip), min0(npts-1, i+ip)
      temp1 = (gamma/(gamma+one))**iabs(j-i)
      temp = temp + temp1
      fmntr(i) = fmntr(i) + rw(j)*temp1
   20 continue
      fmntr(i) = fmntr(i)/temp
   30 continue
c
      return
      end
c
c***********************************************************************
c***********************************************************************
c
      subroutine deffnc(x, k, npde, npts, xmesh, u, ux, uval, uxval)
c
c...this subroutine computes the interpolatory value (uval) and first
c...spatial derivative (uxval) at any point (x) for the k_th component
c...of the solution by using the mesh points (xmesh) and corresponding
c...mesh values for u and ux
c
      implicit real*8 (a-h, o-z)
      dimension xmesh(npts), u(npts,npde), ux(npts,npde)
c
      do 10 i1 = 1, npts-1
      if (xmesh(i1).le.x .and. x.le.xmesh(i1+1)) then
         i = i1
         go to 20
      endif
   10 continue
      if (x.le.xmesh(2)) i = 1
      if (x.ge.xmesh(npts-1)) i = npts-1
   20 h = xmesh(i+1)-xmesh(i) 
      s = (x-xmesh(i))/h
      uval = u(i,k)*fncshp(1,0,s)+ux(i,k)*fncshp(2,0,s)*h
     ++u(i+1,k)*fncshp(3,0,s)+ux(i+1,k)*fncshp(4,0,s)*h
      uxval = (u(i,k)*fncshp(1,1,s)+ux(i,k)*fncshp(2,1,s)*h
     ++u(i+1,k)*fncshp(3,1,s)+ux(i+1,k)*fncshp(4,1,s)*h)/h
      return
      end
c
c***********************************************************************
c***********************************************************************
c
      subroutine slnout(npde, npts, t, y, ydot, touta, ntouta,
     +                  tstep, istop, index, nts, u, ux, uxx, ut,
     +                  uxt, xmesh, xmesht, uu, uux, uut)
c
c...this subroutine outputs the solution values x_t, u, u_x,
c...and u_t at point x and time t in the format prescribed by 
c...subroutine defout
c  
      implicit real*8 (a-h, o-z)
      dimension y((2*npde+1)*npts), ydot((2*npde+1)*npts)
      dimension touta(ntouta), u(npde), ux(npde), uxx(npde)
      dimension ut(npde), uxt(npde), xmesh(npts), xmesht(npts)
      dimension uu(npts, npde), uux(npts, npde), uut(npts, npde)
c
      m = 2*npde+1
      do 20 i = 1, npts
c
c...compute u, ux, ut at time t and x = x_i
c
      x = y(m*i)
      i1 = min0(i, npts-1)
      call drvtvs(npde, npts, i1, x, y, ydot, u, ux, uxx, ut, uxt)
c
c...store x_i, xdot_i, u(k), ux(k), ut(k) in xmesh(i), xmesht(i),
c...uu(i,k), uux(i,k) and uut(i,k), respectively
c
      xmesh(i) = x
      xmesht(i) = ydot(m*i)
      do 10 k = 1, npde
      uu(i,k) = u(k)
      uux(i,k) = ux(k)
      uut(i,k) = ut(k)
   10 continue
   20 continue
c
c...output the solution values and current time stepsize tstep
c 
      call defout(npde, npts, t, xmesh, xmesht, uu, uux, uut, tstep,
     +            touta, ntouta, istop, index, nts) 
c
       return
       end
c
c***********************************************************************
c***********************************************************************
c
      real*8 function zermch(dum)
c
c...this routine computes the machine unit roundoff (zero) 
c...in double precision. it is the smallest positive machine 
c...number such that  1.d0 + zermch .ne. 1.d0
c
      implicit real*8 (a-h, o-z) 
c
      data half/0.5d0/, one/1.d0/
c
      zermch = one 
   10 temp = zermch*half
      if (temp+one.ne.one) then
         zermch = temp
         go to 10
      endif
      return
      end