      subroutine epsode (diffun, pederv, n, t0, h0, y0, tout, eps,
     1                   ierror, mf, index)
c this is the june 24, 1975  version of
c episode..  experimental package for integration of
c systems of ordinary differential equations,
c    dy/dt = f(y,t),  y = (y(1),y(2),...,y(n)) transpose,
c given the initial value of y.
c this code is for the ibm 370/195 at argonne national laboratory
c and is a modification of earlier versions by g.d.byrne
c and a.c.hindmarsh.
c                 references
c 1.  g. d. byrne and a. c. hindmarsh, a polyalgorithm for the
c       numerical solution of ordinary differential equations,
c       ucrl-75652, lawrence livermore laboratory, p. o. box 808,
c       livermore, ca 94550, april 1974. also in acm transactions
c       on mathematical software, 1 (1975), pp. 71-96.
c
c 2.  a. c. hindmarsh and g. d. byrne, episode.. an experimental
c       package for the integration of systems of ordinary
c       differential equations, ucid-30112, l.l.l., may, 1975.
c
c 3.  a. c. hindmarsh, gear.. ordinary differential equation
c       system solver, ucid-30001, rev. 3, l.l.l., december, 1974.
c
c-----------------------------------------------------------------------
c epsode is a driver subroutine for the episode package.
c epsode is to be called once for each output value of t.
c it then makes repeated calls to the core integrator
c subroutine, tstep.
c
c the input parameters are as follows.
c   diffun=  name of double precision subroutine diffun(n,t,y,ydot)
c            which the user must supply and declare external. see
c            note below.
c   pederv=  name of double precision subroutine pederv(n,t,y,pd,n0)
c            which the user must supply and declare external. see
c            note below.
c   n     =  the number of differential equations (used only on
c              first call, unless index = -1).  n must never be
c              increased during a given problem.
c   t0    =  the initial value of t, the independent variable
c              (used for input only on first call).
c   h0    =  the step size in t (used for input only on the
c              first call, unless index = 3 on input).  when
c              index = 3, h0 is the maximum absolute value of
c              the step size to be used.
c   y0    =  a vector of length n containing the initial values of
c              y (used for input only on first call).
c   tout  =  the value of t at which output is desired next.
c              integration will normally go beyond tout and
c              interpolate to t = tout.  (used only for input.)
c   eps   =  the relative error bound (used only on first call,
c              unless index = -1).  this bound is used as follows.
c              let r(i) denote the estimated relative local error
c              in y(i), i.e. the error relative to ymax(i), as
c              measured per step (of size h) or per ss units of t.
c              then eps is a bound on the root-mean-square norm
c              of the vector r, i.e.
c                      n
c              sqrt ( sum ( r(i)**2 )/n ) .lt. eps.
c                     i=1
c              the vector ymax is computed in epsode as described
c              under ierror below.
c              if error control per ss units of t is desired, set ss
c              to a positive number after statement 10 (where it is
c              now set to zero) and update it after statement 60.
c              see also the comments on ss and ymax below.
c  ierror =  the error flag with values and meanings as follow.
c         1    absolute error is controlled.  ymax(i) = 1.0.
c         2    error relative to abs(y) is controlled.  if y(i) = 0.0
c              a divide error will occur.  ymax(i) = abs(y(i)).
c         3    error relative to the largest value of abs(y(i)) seen
c              so far is controlled.  if the initial value of y(i) is
c              0.0, then ymax(i) is set to 1.0 initially and remains
c              at least 1.0.
c   mf    =  the method flag (used only on first call, unless
c              index = -1).  allowed values are 10, 11, 12, 13,
c              20, 21, 22, 23.  mf is an integer with two decimal
c              digits, meth and miter (mf = 10*meth + miter).  (mf
c              can be thought of as the ordered pair (meth,miter).)
c              meth is the basic method indicator.
c                meth = 1 indicates variable-step size, variable-
c                         order adams method, suitable for non-
c                         stiff problems.
c                meth = 2 indicates variable-step size, variable-
c                         order backward differentiation method,
c                         suitable for stiff problems.
c              miter indicates the method of iterative correction
c                (nonlinear system solution).
c                miter = 0 indicates functional iteration (no
c                          partial derivatives needed).
c                miter = 1 indicates a chord or semi-stationary
c                          newton method with closed form (exact)
c                          jacobian, which is computed in the
c                          user supplied subroutine
c                          pederv(n,t,y,pd,n0) described below.
c                miter = 2 indicates a chord  or semi-stationary
c                          newton method with an internally
c                          computed finite difference approximation
c                          to the jacobian.
c                miter = 3 indicates a chord or semi-stationary
c                          newton method with an internally
c                          computed diagonal matrix approximation
c                          to the jacobian, based on a directional
c                          derivative.
c   index =  integer used on input to indicate type of call,
c              with the following values and meanings..
c         1    this is the first call for this problem.
c         0    this is not the first call for this problem,
c              and integration is to continue.
c        -1    this is not the first call for the problem,
c              and the user has reset n, eps, and/or mf.
c         2    same as 0 except that tout is to be hit
c              exactly (no interpolation is done).
c              assumes tout .ge. the current t.
c         3    same as 0 except control returns to calling
c              program after one step.  tout is ignored.
c            since the normal output value of index is 0,
c            it need not be reset for normal continuation.
c
c after the initial call, if a normal return occurred and a normal
c continuation is desired, simply reset tout and call again.
c all other parameters will be ready for the next call.
c a change of parameters with index = -1 can be made after
c either a successful or an unsuccessful return.
c
c the output parameters are as follows.
c   t0    =  the output value of t.  if integration was successful,
c            t0 = tout.  otherwise, t0 is the last value of t
c            reached successfully.
c   h0    =  the step size h used last, whether successfully or not.
c   y0    =  the computed values of y at t = t0.
c   index =  integer used on output to indicate results,
c              with the following values and meanings..
c         0    integration was completed to tout or beyond.
c        -1    the integration was halted after failing to pass the
c              error test even after reducing h by a factor of
c              1.e10 from its initial value.
c        -2    after some initial success, the integration was
c              halted either by repeated error test failures or
c              by a test on eps.  possibly too much accuracy has
c              been requested, or a bad choice of mf was made.
c        -3    the integration was halted after failing to achieve
c              corrector convergence even after reducing h by a
c              factor of 1.e10 from its initial value.
c        -4    immediate halt because of illegal values of input
c              parameters.  see printed message.
c        -5    index was -1 on input, but the desired changes of
c              parameters were not implemented because tout
c              was not beyond t.  interpolation to t = tout was
c              performed as on a normal return.  to continue,
c              simply call again with index = -1 and a new tout.
c        -6    index was 2 on input, but tout was not beyond t.
c              no action was taken.
c
c in addition to epsode, the following subroutines are used by and
c provided in this package:
c   interp(tout,y,n0,y0)  interpolates to give output values at
c                         t = tout  by using data in the y array.
c   tstep(diffun,pederv,y,n0)  is the core integration subroutine, which
c                integrates over a single step and does associated error
c                control.
c   coset  sets coefficients for use in tstep.
c   adjust(y,n0)  adjusts the history array y on reduction of order.
c   pset(diffun,pederv,y,n0,con,miter,ier)  computes and processes the
c                             jacobian matrix, j = df/dy.
c   dec(n,n0,a,ip,ier)  performs the lu decomposition of a matrix.
c   sol(n,n0,a,b,ip)  solves a linear system a*x = b, after dec
c                       has been called for the matrix a.
c note:  pset, dec, and sol are called if and only if miter = 1
c        or miter = 2.
c
c the user must furnish the following double precision subroutines
c      and declare them external in his calling program:
c   diffun(n,t,y,ydot)  computes the function  ydot = f(y,t),
c                       the right hand side of the ordinary
c                       differential equation system, where y
c                       and ydot are vectors of length n.
c   pederv(n,t,y,pd,n0)  computes the n by n jacobian matrix of
c                        partial derivatives and stores it in pd as
c                        an n0 by n0 array.   pd(i,j) is to be set
c                        to the partial derivative of ydot(i) with
c                        respect to y(j).  pederv is called if and
c                        only if miter = 1. for other values of
c                        miter, pederv can be a dummy subroutine.
c
c caution:  at the present time the maximum number of differential
c           equations, which can be solved by episode, is 20.  to
c           change this number to a new value, say nmax, change
c           y(20,13) to y(nmax,13), ymax(20) to ymax(nmax),
c           error(20) to error(nmax), save1(20) to save1(nmax),
c           save2(20) to save2(nmax), pw(400) to pw(nmax*nmax),
c           and ipiv(20) to ipiv(nmax) in the common and dimension
c           statements below.  also change the argument in the
c           if...go to 440 statement (after the common statements)
c           from 20 to nmax.  no other changes need to be made to
c           any other subroutine in this package when the maximum
c           number of equations is changed.  elsewhere, the column
c           length of the y array is n0 instead of 20.  the row
c           length of y can be reduced from 13 to 6 if meth = 2.
c           the array ipiv is used if and only if miter = 1 or
c           miter = 2.  the size of the pw array can be reduced
c           to 1 if miter = 0 or to n if miter = 3.
c
c the common block epcom2 can be accessed externally by the user,
c if he desires to modify the error control behaviour of the package.
c the block contains ymax(20) used to control the error. see above
c and reference (2).
c
c the common block epcom9 can be accessed externally by the user,
c if he desires.  it contains the step size last used successfully
c (hused), the order last used successfully (nqused), the
c number of steps taken so far (nstep), the number of function
c evaluations (diffun calls) so far (nfe), and the number of
c jacobian evaluations so far (nje).
c
c in a data statement below, lout is set to the logical unit number
c for the output of messages during integration.  currently, lout
c = 6.
c-----------------------------------------------------------------------
c
      external diffun, pederv
      integer ierror, index, mf, n
      integer ipiv, jstart, kflag, mfc, nc, nfe, nje,
     1        nqused, nsq, nstep
      integer i, kgo, nhcut, n0
      integer lout
      double precision eps, h0, tout, t0, y0
      double precision epsc, epsj, error, hmax, h, hmin, hused,
     1                 pw, save1, save2, ss, t, uround, ymax
      double precision ayi, d, epsj, h0, t0p, y
      double precision hcut
      double precision four, hundrd, one, ten, zero
      dimension y(20,13)
      dimension y0(n)
c
      common /epcom1/t,h,hmin,hmax,epsc,ss,uround,nc,mfc,kflag,jstart
      common /epcom2/ ymax(20)
      common /epcom3/ error(20)
      common /epcom4/ save1(20)
      common /epcom5/ save2(20)
      common /epcom6/ pw(400)
      common /epcom7/ ipiv(20)
      common /epcom8/ epsj,nsq
      common /epcom9/ hused,nqused,nstep,nfe,nje
c
      data lout /6/
      data hcut /0.1d0/
      data four /4.0d0/, hundrd /1.0d2/, one  /1.0d0/,
     1     ten  /1.0d1/, zero   /0.0d0/
      if (index .eq. 0) go to 20
      if (index .eq. 2) go to 25
      if (index .eq. -1) go to 30
      if (index .eq. 3) go to 40
      if (index .ne. 1) go to 430
      if (eps .le. zero) go to 400
      if (n .le. 0) go to 410
      if (n .gt. 20) go to 440
      if ((t0-tout)*h0 .ge. zero) go to 420
c-----------------------------------------------------------------------
c if initial values for ymax other than those below are desired,
c they should be set here.  all ymax(i) must be positive.  if
c values for hmin or hmax, the bounds on the absolute value of h,
c other than those below, are desired, they also should be set here.
c if error per ss units of t is to be controlled, ss should be set
c to a positive value below.  error per unit step is controlled
c when ss = 1.  the default value for ss is 0 and yields control
c of error per step.
c-----------------------------------------------------------------------
c set uround, the machine roundoff constant, here.
c use statement below for short precision on ibm 360 or 370.
c     uround = 9.53674e-7
c use statement below for single precision on cdc 7600 or 6600.
c     uround = 7.105427406e-15
c use statement below for long precision on ibm 360 or 370.
c     uround = 2.220446047d-16
      uround = d1mach(4)
      do 10 i = 1,n
        go to (5, 6, 7), ierror
c ierror   =   1, 2, 3  ------------------------------------------------
  5     ymax(i) = one
        go to 10
  6     ymax(i) = dabs(y0(i))
        go to 10
  7     ymax(i) = dabs(y0(i))
        if (ymax(i) .eq. zero) ymax(i) = one
 10     y(i,1) = y0(i)
      nc = n
      t = t0
      h = h0
      if ((t+h) .eq. t) write(lout,15) t
 15   format(/47h---  message from subroutine epsode in episode,,
     1       24h the o.d.e. solver.  ---/22h warning.. t + h = t =,
     2       e18.8,18h in the next step./)
      hmin = dabs(h0)
      hmax = dabs(t0 - tout)*ten
      epsc = eps
      mfc = mf
      jstart = 0
      ss = zero
      n0 = n
      nsq = n0*n0
      epsj = dsqrt(uround)
      nhcut = 0
      go to 50
c t0p is the previous output value of t0 for use in hmax. --------------
  20  hmax = dabs(tout - t0p)*ten
      go to 80
  25  hmax = dabs(tout - t0p)*ten
c
       if ((t-tout)*h .ge. zero) go to 460
       go to 85
c
 30   if ((t-tout)*h .ge. zero) go to 450
      if (mf .ne. mfc) jstart = -1
      nc = n
      epsc = eps
      mfc = mf
      go to 45
c
 40   hmax = h0
c
 45   if ((t+h) .eq. t) write(lout,15) t
c
 50   call tstep (diffun, pederv, y, n0)
c
      kgo = 1 - kflag
      go to (60, 100, 200, 300), kgo
c kflag  =   0,  -1,  -2,  -3  -----------------------------------------
c
 60   continue
c-----------------------------------------------------------------------
c normal return from tstep.
c
c the weights ymax(i) are updated.  if different values are desired,
c they should be set here.  if ss is to be updated for control of
c error per ss units of t, it should also be done here.  a test is
c made to determine if eps is too small for machine precision.
c
c any other tests or calculations that are required after each step
c should be inserted here.
c
c if index = 3, y0 is set to the current y values on return.
c if index = 2, h is controlled to hit tout (within roundoff
c error), and then the current y values are put in y0 on
c return.  for any other value of index, control returns to
c the integrator unless tout has been reached.  then
c interpolated values of y are computed and stored in y0 on
c return.
c if interpolation is not desired, the call to interp should
c be deleted and control transferred to statement 500 instead
c of 520.
c-----------------------------------------------------------------------
      d = zero
      do 70 i = 1,n
        ayi = dabs(y(i,1))
        go to (70, 66, 67), ierror
c ierror  =     1,  2,  3  ---------------------------------------------
 66     ymax(i) = ayi
        go to 70
  67    ymax(i) = dmax1(ymax(i), ayi)
  70    d = d + (ayi/ymax(i))**2
      d = d*(uround/eps)**2
      if (d .gt. dfloat(n)) go to 250
      if (index .eq. 3) go to 500
      if (index .eq. 2) go to 85
 80   if ((t-tout)*h .lt. zero) go to 45
      call interp (tout, y, n0, y0)
      t0 = tout
      go to 520
 85   if (((t+h)-tout)*h .le. zero) go to 45
      if (dabs(t-tout) .le. hundrd*uround*hmax) go to 500
      if ((t-tout)*h .ge. zero) go to 500
      h = (tout - t)*(one - four*uround)
      jstart = -1
      go to 45
c-----------------------------------------------------------------------
c on an error return from tstep, an immediate return occurs if
c kflag = -2, and recovery attempts are made otherwise.
c to recover, h and hmin are reduced by a factor of .1 up to 10
c times before giving up.
c-----------------------------------------------------------------------
 100  write (lout,101)
 101  format(/47h---  message from subroutine epsode in episode,,
     1        24h the o.d.e. solver.  ---/)
      write(lout,105) t,hmin
 105  format(//35h kflag = -1 from integrator at t = ,e18.8/
     1       40h  error test failed with abs(h) = hmin =,e18.8/)
 110  if (nhcut .eq. 10) go to 150
      nhcut = nhcut + 1
      hmin = hcut*hmin
      h = hcut*h
      write (lout,115) h
 115  format(24h  h has been reduced to ,e18.8,
     1       26h  and step will be retried//)
      jstart = -1
      go to 45
c
 150  write (lout,155)
 155  format(//44h problem appears unsolvable with given input//)
      go to 500
c
 200  write (lout,101)
      write (lout,205) t,h,eps
 205  format(//16h kflag = -2  t =,e18.8,4h h =,e18.8,6h eps =,e18.8/
     1       50h  the requested error is too small for integrator.//)
      go to 500
c
 250  write (lout,101)
      write (lout,255) t,eps
 255  format(//47h integration halted by subroutine epsode at t =,
     1       e18.8/43h eps is too small for machine precision and/
     2       29h problem being solved.  eps =,e18.8//)
      kflag = -2
      go to 500
c
 300  write (lout,101)
      write (lout,305) t
 305  format(//34h kflag = -3 from integrator at t =,e18.8/
     1       45h  corrector convergence could not be achieved/)
      go to 110
c
 400  write (lout,101)
      write (lout,405) eps
 405  format(//35h illegal input.. eps .le. 0. eps = ,e18.8//)
      index = -4
      return
c
 410  write (lout,101)
      write (lout,415) n
 415  format(//31h illegal input.. n .le. 0. n = ,i8//)
      index = -4
      return
c
 420  write (lout,101)
      write (lout,425) t0,tout,h0
 425  format(//39h illegal input.. (t0 - tout)*h0 .ge. 0./
     1       5h t0 =,e18.8,7h tout =,e18.8,5h h0 =,e18.8//)
      index = -4
      return
c
 430  write (lout,101)
      write (lout,435) index
 435  format(//24h illegal input.. index =,i8//)
      index = -4
      return
c
 440  write (lout,101)
      write (lout,445) n
 445  format (//39h illegal input.  the number of ordinary/
     1        43h differential equations being solved is n =, i6/
     2        43h storage allocation in subroutine epsode is/
     3        47h too small.  see writeups#epsode public......../)
      index = -4
      return
c
 450  write (lout,101)
      write (lout,455) t,tout,h
 455  format(//46h index = -1 on input with (t - tout)*h .ge. 0./
     1       44h interpolation was done as on normal return./
     2       41h desired parameter changes were not made./
     3       4h t =,e18.8,7h tout =,e18.8,4h h =,e18.8//)
      call interp (tout, y, n0, y0)
      t0 = tout
      index = -5
      return
c
 460  write (lout,101)
      write (lout,465) t,tout,h
 465  format(//45h index = 2 on input with (t - tout)*h .ge. 0./
     1       4h t =,e18.8,7h tout =,e18.8,4h h =,e18.8//)
      index = -6
      return
c
 500  t0 = t
      do 510 i = 1,n
 510    y0(i) = y(i,1)
 520  index = kflag
      t0p = t0
      h0 = hused
      if (kflag .ne. 0) h0 = h
      return
c-----------------------  end of subroutine epsode ---------------------
      end
      subroutine interp (tout, y, n0, y0)
c-----------------------------------------------------------------------
c subroutine interp computes interpolated values of the dependent
c variable y and stores them in y0.  the interpolation is to the
c point t = tout and uses the nordsieck history array y as follows..
c                             nq
c                  y0(i)  =  sum  y(i,j+1)*s**j ,
c                            j=0
c where s = -(t-tout)/h.
c-----------------------------------------------------------------------
c caution:  not all members of epcom1 are used in this subroutine.
c-----------------------------------------------------------------------
      integer n0
      integer jstart, kflag, mf, n
      integer i, j, l
      double precision tout, y, y0
      double precision eps, h, hmax, hmin, ss, t, uround
      double precision s, s1
      double precision one
      dimension y0(n0),y(n0,13)
c
      common /epcom1/ t,h,hmin,hmax,eps,ss,uround,n,mf,kflag,jstart
      data one /1.0d0/
      do 10 i = 1,n
 10     y0(i) = y(i,1)
      l = jstart + 1
      s = (tout - t)/h
      s1 = one
      do 30 j = 2,l
        s1 = s1*s
        do 20 i = 1,n
 20       y0(i) = y0(i) + s1*y(i,j)
 30     continue
      return
c-----------------------  end of subroutine interp  --------------------
      end
      subroutine tstep (diffun, pederv, y, n0)
c-----------------------------------------------------------------------
c tstep performs one step of the integration of an initial value
c problem for a system of ordinary differential equations.
c communication with tstep is via the following variables..
c
c   y       an n0 by lmax array containing the dependent variables
c             and their scaled derivatives.  lmax is currently 6 for
c             the variable step backward differentiation formulas,
c             and 13 for the variable step adams formulas.
c             (lmax -1) = maxder, the maximum order used.
c             see subroutine coset.  y(i,j+1) contains the
c             j-th derivative of y(i), scaled by h**j/factorial(j)
c             for j = 0,1,...,nq, where nq is the current order.
c   n0      a constant integer .ge. n, used for dimensioning
c             purposes.
c   t       the independent variable, updated on each step taken.
c   h       the step size to be attempted on the next step.
c             h is altered by the error control algorithm during
c             the solution of the problem. h can be either positive
c             or negative, but its sign must remain constant
c             throughout the problem run.
c   hmin,   the minimum and maximum absolute values of the step
c    hmax     size to be used for the step. these may be changed at
c             any time, but the change will not take effect until the
c             next change in h is made.
c   eps     the relative error bound.  see description in
c             subroutine epsode.
c   ss      the size of the time interval to be used for error
c             control.  a default value of 0 is used to produce
c             control of error per step.  see subroutine epsode.
c   uround  the unit of roundoff for the computer being used.
c   n       the number of first order ordinary differential
c             equations being solved.
c   mf      the method flag.  see description in subroutine epsode.
c   kflag   a completion code with the following meanings..
c                     0  the step was succesful.
c                    -1  the requested error could not be achieved
c                          with abs(h) = hmin.
c                    -2  the requested error is smaller than can
c                          be handled for this problem.
c                    -3  corrector convergence could not be
c                          achieved for abs(h) = hmin.
c             on a return with kflag negative, the values of t and
c             the y array are as of the beginning of the last
c             step and h is the last step size attempted.
c   jstart  an integer used on input and output.
c             on input, it has the following values and meanings..
c                     0  perform the first step.
c                 .gt.0  take a new step continuing from the last.
c                 .lt.0  take the next step with a new value of
c                          h and/or mf.
c             on exit, jstart is set to nq, the current order of the
c             method.
c   ymax      an array of n elements with which the estimated local
c               errors in y are compared.
c   error     an array of n elements.  error(i)/tq(2) is the
c               estimated local error in y(i) per ss units of
c               t or per step (of size h).
c   save1,    two arrays for working storage,
c     save2     each of length n.
c   pw        a block of locations used for the partial derivatives
c               of f with respect to y, if miter is not o.  see
c               description in subroutine epsode.
c   ipiv      an integer array of length n, which is used for pivot
c               information for the linear algebraic system in the
c               correction process, when miter = 1 or 2.
c
c the common block epcm10, declared below, is primarily intended
c for internal use, but it can be accessed externally.
c-----------------------------------------------------------------------
      external pederv
      integer n0
      integer ipiv, jstart, kflag, l, lmax, meth, mf, n, nfe, nje,
     1        nq, nqindx, nqused, nstep
      integer i, iback, ier, iredo, j, j1, j2, m, mfold, mio,
     1        miter, miter1, newj, nstepj
      integer istepj, kfc, kfh, maxcor
      double precision y
      double precision el, eps, error, h, hmax, hmin, hused, pw,
     1                 save1, save2, ss, t, tau, tq, uround, ymax
      double precision bnd, cnquot, con, conp, crate, d, drc,
     1                 d1, e, edn, eta, etamax, etamin, etaq, etaqm1,
     2                 etaqp1, eup, flotl, flotn, hold, hrl1, phrl1,
     3                 prl1, r, rc, rl1, r0, r1, told
      double precision addon, bias1, bias2, bias3, crdown, delrc,
     1                 etacf, etamin, etamxf, etamx1, etamx2,
     2                 etamx3, onepsm, short, thresh
      double precision one, pt5, zero
      dimension y(n0,13)
c
      common /epcom1/ t,h,hmin,hmax,eps,ss,uround,n,mf,kflag,jstart
      common /epcom2/ ymax(1)
      common /epcom3/ error(1)
      common /epcom4/ save1(1)
      common /epcom5/ save2(1)
      common /epcom6/ pw(1)
      common /epcom7/ ipiv(1)
      common /epcom9/ hused,nqused,nstep,nfe,nje
      common /epcm10/ tau(13),el(13),tq(5),lmax,meth,nq,l,nqindx
c
      data istepj /20/, kfc /-3/, kfh /-7/, maxcor /3/
      data addon  /1.0d-6/,    bias1  /2.5d1/, bias2  /2.5d1/,
     1     bias3  /1.0d2/,     crdown /0.1d0/, delrc  /0.3d0/,
     2     etacf  /0.25d0/,    etamin /0.1d0/, etamxf /0.2d0/,
     3     etamx1 /1.0d4/,     etamx2 /1.0d1/, etamx3 /1.5d0/,
     4     onepsm /1.00001d0/, short  /0.1d0/, thresh /1.3d0/
      data one /1.0d0/, pt5 /0.5d0/, zero /0.0d0/
      kflag = 0
      told = t
      flotn = dfloat(n)
      if (jstart .gt. 0) go to 200
      if (jstart .ne. 0) go to 150
c-----------------------------------------------------------------------
c on the first call, the order is set to 1 and the initial
c derivatives are calculated.  etamax is the maximum ratio by
c which h can be increased in a single step.  it is 1.e04 for the
c first step to compensate for the small initial h, then 10 for
c the next 10 steps, and then 1.5 thereafter.  if a failure
c occurs (in corrector convergence or error test), etamax is set at 1
c for the next increase.  etamin = .1 is the minimum ratio by which
c h can be reduced on any retry of a step.
c-----------------------------------------------------------------------
      call diffun (n, t, y, save1)
      do 110 i = 1,n
 110    y(i,2) = h*save1(i)
      meth = mf/10
      miter = mf - 10*meth
      miter1 = miter + 1
      mfold = mf
      nq = 1
      l = 2
      tau(1) = h
      prl1 = one
      rc = zero
      etamax = etamx1
      nqindx = 2
      nstep = 0
      nstepj = 0
      nfe = 1
      nje = 0
      go to 200
c-----------------------------------------------------------------------
c if the user has changed h, then y must be rescaled.  if the
c user has changed miter, then newj is set to miter to force
c the partial derivativees to be updated, if they are being used.
c-----------------------------------------------------------------------
 150  if (mf .eq. mfold) go to 170
      mio = miter
      meth = mf/10
      miter = mf - 10*meth
      mfold = mf
      if (miter .eq. mio) go to 170
      newj = miter
      miter1 = miter + 1
 170  if (h .eq. hold) go to 200
      eta = h/hold
      h = hold
      iredo = 3
      go to 185
  180 eta = dmax1(eta,hmin/dabs(h),etamin)
  185 eta = dmin1(eta,hmax/dabs(h),etamax)
      r1 = one
      do 190 j = 2,l
        r1 = r1*eta
        do 190 i = 1,n
 190      y(i,j) = y(i,j)*r1
      h = h*eta
      rc = rc*eta
      if (iredo .eq. 0) go to 690
c-----------------------------------------------------------------------
c this section computes the predicted values by effectively
c multiplying the y array by the pascal triangle matrix.  then
c coset is called to obtain el, the vector of coefficients of
c length nq + 1.  rc is the ratio of new to old values of the
c coefficient h/el(2).  when rc differs from 1 by more than
c delrc, newj is set to miter to force the partial derivatives
c to be updated, if used.  delrc is 0.3.  in any case, the partial
c derivatives are updated at least every 20-th step.
c-----------------------------------------------------------------------
 200  t = t + h
      do 210 j1 = 1,nq
        do 210 j2 = j1,nq
          j = (nq + j1) - j2
          do 210 i = 1,n
 210        y(i,j) = y(i,j) + y(i,j+1)
      call coset
      bnd = flotn*(tq(4)*eps)**2
      rl1 = one/el(2)
      rc = rc*(rl1/prl1)
      prl1 = rl1
      if (nstep .ge. nstepj+istepj) newj = miter
      drc = dabs(rc-one)
      if (drc .le. delrc) go to 215
      newj = miter
      crate = one
      rc = one
      go to 220
  215 if ((miter .ne. 0) .and. (drc .ne. zero)) crate = one
c-----------------------------------------------------------------------
c up to 3 corrector iterations are taken.  a convergence test is made
c on the root mean square norm of each correction, using bnd, which
c is dependent on eps.  the sum of the corrections is accumulated in
c the vector error.  the y array is not altered in the corrector
c loop.  the updated y vector is stored temporarily in save1.
c-----------------------------------------------------------------------
 220  do 230 i = 1,n
 230    error(i) = zero
      m = 0
      call diffun (n, t, y, save2)
      nfe = nfe + 1
      if (newj .le. 0) go to 290
c-----------------------------------------------------------------------
c if indicated, the matrix p = i - h*rl1*j is reevaluated before
c starting the corrector iteration.  newj is set to 0 as an
c indicator that this has been done.  if miter = 1 or 2, p is
c computed and processed in pset.  if miter = 3, the matrix  is
c p = i - h*rl1*d, where d is a diagonal matrix.  rl1 is 1/el(2).
c-----------------------------------------------------------------------
      newj = 0
      rc = one
      nje = nje + 1
      nstepj = nstep
      go to (250, 240, 260), miter
 240  nfe = nfe + n
 250  con = -h*rl1
      call pset(diffun, pederv, y, n0, con, miter, ier)
      if (ier .ne. 0) go to 420
      go to 350
 260  r = rl1*short
      do 270 i = 1,n
 270    pw(i) = y(i,1) + r*(h*save2(i) - y(i,2))
      call diffun(n, t, pw, save1)
      nfe = nfe + 1
      hrl1 = h*rl1
      do 280 i = 1,n
        r0 = h*save2(i) - y(i,2)
        pw(i) = one
        d = short*r0 - h*(save1(i) - save2(i))
        save1(i) = zero
        if (dabs(r0) .lt. uround*ymax(i)) go to 280
        if (dabs(d) .eq. zero) go to 420
        pw(i) = short*r0/d
        save1(i) = pw(i)*rl1*r0
 280    continue
      go to 370
 290  go to (295, 350, 350, 310), miter1
c-----------------------------------------------------------------------
c in the case of functional iteration, y is updated directly from
c the result of the last diffun call.
c-----------------------------------------------------------------------
 295  d = zero
      do 300 i = 1,n
        r = rl1*(h*save2(i) - y(i,2))
        d = d + ((r - error(i))/ymax(i))**2
        save1(i) = y(i,1) + r
 300    error(i) = r
      go to 400
c-----------------------------------------------------------------------
c in the case of a chord method, the residual -g(y sub n(m))
c is computed and the linear system with that as right-hand side
c and p as coefficient matrix is solved, using the lu decomposition
c of p if miter = 1 or 2.  if miter = 3 the scalar h*rl1 is updated.
c-----------------------------------------------------------------------
 310  phrl1 = hrl1
      hrl1 = h*rl1
      if (hrl1 .eq. phrl1) go to 330
      r = hrl1/phrl1
      do 320 i = 1,n
        d = one - r*(one - one/pw(i))
        if (dabs(d) .eq. zero) go to 440
 320    pw(i) = one/d
 330  do 340 i = 1,n
 340    save1(i) = pw(i)*(rl1*h*save2(i) - (rl1*y(i,2) + error(i)))
      go to 370
 350  do 360 i = 1,n
 360    save1(i) = rl1*h*save2(i) - (rl1*y(i,2) + error(i))
      call sol (n, n0, pw, save1, ipiv)
 370  d = zero
      do 380 i = 1,n
        error(i) = error(i) + save1(i)
        d = d + (save1(i)/ymax(i))**2
 380    save1(i) = y(i,1) + error(i)
c-----------------------------------------------------------------------
c test for convergence.  if m .gt. 0, an estimate of the square of
c the convergence rate constant is stored in crate, and this is used
c in the test.
c-----------------------------------------------------------------------
  400 if (m .ne. 0) crate = dmax1(crdown*crate,d/d1)
      if (d*dmin1(one,crate) .le. bnd) go to 450
      d1 = d
      m = m + 1
      if (m .eq. maxcor) go to 410
      call diffun (n, t, save1, save2)
      go to (295, 350, 350, 310), miter1
c-----------------------------------------------------------------------
c the corrector iteration failed to converge in 3 tries. if partial
c derivatives are involved but are not up to date, they are
c reevaluated for the next try.  otherwise the y array is restored
c to its values before prediction, and h is reduced,
c if possible.  if not, a no-convergence exit is taken.
c-----------------------------------------------------------------------
 410  nfe = nfe + maxcor - 1
      if (newj .eq. -1) go to 440
 420  t = told
      etamax = one
      do 430 j1 = 1,nq
        do 430 j2 = j1,nq
          j = (nq + j1) - j2
          do 430 i = 1,n
 430        y(i,j) = y(i,j) - y(i,j+1)
      if (dabs(h) .le. hmin*onepsm) go to 680
      eta = etacf
      iredo = 1
      go to 180
 440  newj = miter
      go to 220
c-----------------------------------------------------------------------
c the corrector has converged.  newj is set to -1 if partial
c derivatives were used, to signal that they may need updating on
c subsequent steps.  the error test is made and control passes to
c statement 500 if it fails.
c-----------------------------------------------------------------------
 450  if (miter .ne. 0) newj = -1
      nfe = nfe + m
      d = zero
      do 460 i = 1,n
 460    d = d + (error(i)/ymax(i))**2
      e = flotn*(tq(2)*eps)**2
      if (d .gt. e) go to 500
c-----------------------------------------------------------------------
c after a successful step, the y array, tau, nstep, and nqindx are
c updated, and a new value of h at order nq is computed.
c the vector tau contains the nq + 1 most recent values of h.
c a change in nq up or down by 1 is considered if nqindx = 0.
c if nqindx = 1 and nq .lt. maxder, then error is saved
c for use in a possible order increase on the next step.
c a change in h or nq is made only of the increase in h
c is by a factor of at least 1.3.
c if not, nqindx is set to 2 to prevent testing for that many
c steps.  if nq is changed, nqindx is set to nq + 1 (new value).
c-----------------------------------------------------------------------
      kflag = 0
      iredo = 0
      nstep = nstep + 1
      hused = h
      nqused = nq
      do 470 iback = 1,nq
        i = l - iback
 470    tau(i+1) = tau(i)
      tau(1) = h
      do 480 j = 1,l
        do 480 i = 1,n
 480      y(i,j) = y(i,j) + error(i)*el(j)
      nqindx = nqindx - 1
      if ((l .eq. lmax) .or. (nqindx .ne. 1)) go to 495
      do 490 i = 1,n
 490    y(i,lmax) = error(i)
      conp = tq(5)
 495  if (etamax .ne. one) go to 520
      if (nqindx .lt. 2) nqindx = 2
      go to 690
c-----------------------------------------------------------------------
c the error test failed.  kflag keeps track of multiple failures.
c t and the y array are restored to their previous values.  a new
c h for a retry of the step is computed.  the order is kept fixed.
c-----------------------------------------------------------------------
 500  kflag = kflag - 1
      t = told
      do 510 j1 = 1,nq
        do 510 j2 = j1,nq
          j = (nq + j1) - j2
          do 510 i = 1,n
 510        y(i,j) = y(i,j) - y(i,j+1)
      newj = miter
      etamax = one
      if (dabs(h) .le. hmin*onepsm) go to 660
      if (kflag .le. kfc) go to 630
      iredo = 2
c compute ratio of new h to current h at the current order. ------------
  520 flotl = dfloat(l)
      etaq = one/((bias2*d/e)**(pt5/flotl) + addon)
      if ((nqindx .ne. 0) .or. (kflag .ne. 0)) go to 580
      etaqm1 = zero
      if (nq .eq. 1) go to 540
c compute ratio of new h to current h at the current order less one. ---
      d = zero
      do 530 i = 1,n
 530    d = d + (y(i,l)/ymax(i))**2
      edn = flotn*(tq(1)*eps)**2
      etaqm1 = one/((bias1*d/edn)**(pt5/(flotl - one)) + addon)
 540  etaqp1 = zero
      if (l .eq. lmax) go to 560
c compute ratio of new h to current h at current order plus one. -------
      cnquot = (tq(5)/conp)*(h/tau(2))**l
      d = zero
      do 550 i = 1,n
 550    d = d + ((error(i) - cnquot*y(i,lmax))/ymax(i))**2
      eup = flotn*(tq(3)*eps)**2
      etaqp1 = one/((bias3*d/eup)**(pt5/(flotl + one)) + addon)
 560  nqindx = 2
      if (etaq .ge. etaqp1) go to 570
      if (etaqp1 .gt. etaqm1) go to 600
      go to 590
 570  if (etaq .lt. etaqm1) go to 590
 580  if ((etaq .lt. thresh) .and. (kflag .eq. 0)) go to 690
      eta = etaq
      if ((kflag .le. -2) .and. (eta .gt. etamxf)) eta = etamxf
      go to 180
 590  if (etaqm1 .lt. thresh) go to 690
      call adjust (y, n0)
      l = nq
      nq = nq - 1
      eta = etaqm1
      nqindx = l
      go to 180
 600  if (etaqp1 .lt. thresh) go to 690
      nq = l
      eta = etaqp1
      l = l + 1
      do 610 i = 1,n
 610    y(i,l) = zero
      nqindx = l
      go to 180
c-----------------------------------------------------------------------
c control reaches this section if 3 or more consecutive failures
c have occurred.  it is assumed that the elements of the y array
c have accumulated errors of the wrong order.  the order is reduced
c by one, if possible.  then h is reduced by a factor of 0.1 and
c the step is retried.  after a total of 7 consecutive failures,
c an exit is taken with kflag = -2.
c-----------------------------------------------------------------------
 630  if (kflag .eq. kfh) go to 670
      if (nq .eq. 1) go to 640
      eta = etamin
      call adjust (y, n0)
      l = nq
      nq = nq - 1
      nqindx = l
      go to 180
  640 eta = dmax1(etamin,hmin/dabs(h))
      h = h*eta
      call diffun (n, t, y, save1)
      nfe = nfe + 1
      do 650 i = 1,n
 650    y(i,2) = h*save1(i)
      nqindx = 10
      go to 200
c-----------------------------------------------------------------------
c all returns are made through this section.  h is saved in hold
c to allow the caller to change h on the next step.
c-----------------------------------------------------------------------
 660  kflag = -1
      go to 700
 670  kflag = -2
      go to 700
 680  kflag = -3
      go to 700
 690  etamax = etamx3
      if (nstep .le. 10) etamax = etamx2
 700  hold = h
      jstart = nq
      return
c-----------------------  end of subroutine tstep  ---------------------
      end
      subroutine coset
c-----------------------------------------------------------------------
c coset is called by tstep and sets coefficients for use there.
c
c for each order nq, the coefficients in el are calculated by use of
c  the generating polynomial lambda(x), with coefficients el(i):
c      lambda(x) = el(1) + el(2)*x + ... + el(nq+1)*(x**nq).
c for the backward differentiation formulas,
c                    nq
c      lambda(x) = product (1 + x/xi(i) ) .
c                   i = 1
c for the adams formulas,
c                              nq-1
c      (d/dx) lambda(x) = c * product (1 + x/xi(i) ) ,
c                              i = 1
c      lambda(-1) = 0,    lambda(0) = 1,
c where c is a normalization constant.
c in both cases, xi(i) is defined by
c      h*xi(i) = t sub n  -  t sub (n-i)
c              = h + tau(1) + tau(2) + ... tau(i-1).
c
c coset also sets maxder, the maximum order of the formulas
c available. currently this is 5 for the backward differentiation
c formulas, and 12 for the adams formulas.  to use different
c values (.le. 13),  change the numbers in statements 1 and 2 below.
c
c in addition to variables described previously, communication
c with coset uses the following..
c   tau    = a vector of length 13 containing the past nq values
c            of h.
c   el     = a vector of length 13 in which coset stores the
c            coefficients for the corrector formula.
c   tq     = a vector of length 5 in which coset stores constants
c            used for the convergence test, the error test, and
c            selection of h at a new order.
c   lmax   = maxder + 1, where maxder is the maximum order
c            available.  lmax is the maximum number of columns
c            of the y array to be used.
c   meth   = the basic method indicator.
c   nq     = the current order.
c   l      = nq + 1, the length of the vector stored in el, and
c            the number of columns of the y array being used.
c   nqindx = a counter controlling the frequency of order changes.
c            an order change is about to be considered if
c            nqindx = 1.
c-----------------------------------------------------------------------
c caution:  not all members of epcom1 are used in this subroutine.
c-----------------------------------------------------------------------
      integer jstart, kflag, l, lmax, meth, mf, n, nq, nqindx
      integer i, iback, j, jp1, jstart, kflag, l, lmax, maxder,
     1        meth, mf,n, nq, nqindx, nqm1
      double precision el, eps, h, hmax, hmin, ss, t, tau, tq,
     1                 uround
      double precision ahdss, cnqm1, csum, elp, em, em0, floti,
     1                 flotl, flotnq, hsum, hsum1, prod, rxi, s, xi
      double precision cortes
      double precision one, six, two, zero
      dimension em(13)
c
      common /epcom1/ t,h,hmin,hmax,eps,ss,uround,n,mf,kflag,jstart
      common /epcm10/ tau(13),el(13),tq(5),lmax,meth,nq,l,nqindx
      data cortes /0.1d0/
      data one  /1.0d0/, six /6.0d0/, two /2.0d0/, zero /0.0d0/
      ahdss = one
      if (ss .ne. zero) ahdss = dabs(h)/ss
      flotl = dfloat(l)
      nqm1 = nq - 1
      go to (1, 2), meth
 1    maxder = 12
      go to 100
c
 2    maxder = 5
      go to 200
c
 100  if (nq .ne. 1) go to 110
      el(1) = one
      el(2) = one
      tq(1) = one
      tq(2) = two*ahdss
      tq(3) = six*tq(2)
      tq(5) = one
      go to 300
 110  hsum = h
      em(1) = one
      flotnq = flotl - one
      do 115 i = 2,l
 115    em(i) = zero
      do 150 j = 1,nqm1
        if ((j .ne. nqm1) .or. (nqindx .ne. 1)) go to 130
        s = one
        csum = zero
        do 120 i = 1,nqm1
          csum = csum + s*em(i)/dfloat(i+1)
 120      s = -s
        tq(1) = ahdss*em(nqm1)/(flotnq*csum)
 130    rxi = h/hsum
        do 140 iback = 1,j
          i = (j + 2) - iback
 140      em(i) = em(i) + em(i-1)*rxi
 150    hsum = hsum + tau(j)
c compute integral from -1 to 0 of polynomial and of x times it. -------
      s = one
      em0 = zero
      csum = zero
      do 160 i = 1,nq
        floti = dfloat(i)
        em0 = em0 + s*em(i)/floti
        csum = csum + s*em(i)/(floti+1)
 160    s = -s
c in el, form coefficients of normalized integrated polynomial. --------
      s = one/em0
      el(1) = one
      do 170 i = 1,nq
  170   el(i+1) = s*em(i)/dfloat(i)
      xi = hsum/h
      tq(2) = ahdss*xi*em0/csum
      tq(5) = xi/el(l)
      if (nqindx .ne. 1) go to 300
c for higher order control constant, multiply polynomial by 1+x/xi(q). -
      rxi = one/xi
      do 180 iback = 1,nq
        i = (l + 1) - iback
 180    em(i) = em(i) + em(i-1)*rxi
c compute integral of polynomial. --------------------------------------
      s = one
      csum = zero
      do 190 i = 1,l
      csum = csum + s*em(i)/dfloat(i+1)
 190    s = -s
      tq(3) = ahdss*flotl*em0/csum
      go to 300
c
 200  do 210 i = 3,l
 210    el(i) = zero
      el(1) = one
      el(2) = one
      hsum = h
      hsum1 = zero
      prod = one
      rxi = one
      if (nq .eq. 1) go to 240
      do 230 j = 1,nqm1
c in el, construct coefficients of (1+x/xi(1))*...*(1+x/xi(j+1)). ------
        hsum = hsum + tau(j)
        hsum1 = hsum1 + tau(j)
        prod = prod*(hsum/hsum1)
        rxi = h/hsum
        jp1 = j + 1
        do 220 iback = 1,jp1
          i = (j + 3) - iback
 220      el(i) = el(i) + el(i-1)*rxi
 230    continue
 240  tq(2) = ahdss*el(2)*(one + prod)
      tq(5) = (one + prod)/el(l)
      if (nqindx .ne. 1) go to 300
      cnqm1 = rxi/el(l)
      elp = el(2) - rxi
      tq(1) = ahdss*elp/cnqm1
      hsum = hsum + tau(nq)
      rxi = h/hsum
      elp = el(2) + rxi
      tq(3) = ahdss*elp*rxi*(one + prod)*(flotl + one)
 300  tq(4) = cortes*tq(2)
      lmax = maxder + 1
      return
c-----------------------  end of subroutine coset  ---------------------
      end
      subroutine adjust (y, n0)
c-----------------------------------------------------------------------
c this subroutine adjusts the y array on reduction of order.
c see reference 1 for details.
c-----------------------------------------------------------------------
c caution:  not all members of epcom1 are used in this subroutine.
c-----------------------------------------------------------------------
      integer n0
      integer jstart, kflag, l, lmax, meth, mf, n, nq, nqindx
      integer i, iback, j, jp1, nqm1, nqm2
      double precision y
      double precision el, eps, h, hmax, hmin, ss, t, tau, tq, uround
      double precision hsum, xi
      double precision one, zero
      dimension y(n0,13)
c
      common /epcom1/ t,h,hmin,hmax,eps,ss,uround,n,mf,kflag,jstart
      common /epcm10/ tau(13),el(13),tq(5),lmax,meth,nq,l,nqindx
      data one /1.0d0/, zero /0.0d0/
      if (nq .eq. 2) return
      nqm1 = nq - 1
      nqm2 = nq - 2
      go to (100, 200), meth
c
 100  do 110 j = 1,lmax
 110    el(j) = zero
      el(2) = one
      hsum = zero
      do 130 j = 1,nqm2
c construct coefficients of x*(x+xi(1))*...*(x+xi(j)). -----------------
        hsum = hsum + tau(j)
        xi = hsum/h
        jp1 = j + 1
        do 120 iback = 1,jp1
          i = (j + 3) - iback
 120      el(i) = el(i)*xi + el(i-1)
 130    continue
c construct coefficients of integrated polynomial. ---------------------
      do 140 j = 2,nqm1
  140   el(j+1) = dfloat(nq)*el(j)/dfloat(j)
      go to 300
c
 200  do 210 j = 1,lmax
 210    el(j) = zero
      el(3) = one
      hsum = zero
      do 230 j = 1,nqm2
c construct coefficients of x*x*(x+xi(1))*...*(x+xi(j)). ---------------
        hsum = hsum + tau(j)
        xi = hsum/h
        jp1 = j + 1
        do 220 iback = 1,jp1
          i = (j + 4) - iback
 220      el(i) = el(i)*xi + el(i-1)
 230    continue
c
c subtract correction terms from y array. ------------------------------
 300  do 320 j = 3,nq
        do 310 i = 1,n
 310      y(i,j) = y(i,j) - y(i,l)*el(j)
 320    continue
      return
c-----------------------  end of subroutine adjust  --------------------
      end
      subroutine pset (diffun, pederv, y, n0, con, miter, ier)
c-----------------------------------------------------------------------
c pset is called by tstep to compute and to process the matrix
c p = i - (h/el(2))*j, where j is an approximation to the
c jacobian.  j is computed by either the user supplied
c subroutine pederv, when miter = 1, or by finite differences,
c when miter = 2.  j is stored in pw and replaced by p, using
c con = -h/el(2).  then p is subjected to an lu decomposition
c for later solution of linear algebraic systems with p as the
c coefficient matrix.
c
c in addition to variables described previously, communication
c with pset uses the following..
c   epsj = sqrt(uround), used in the numerical jacobian increments.
c   nsq  = n0**2.
c-----------------------------------------------------------------------
c caution:  not all epcom1 variables are used inthis subroutine.
c-----------------------------------------------------------------------
      integer ier, miter, n0
      integer ipiv, jstart, kflag, mf, n, nsq
      integer i, ier, j, j1, n
      double precision con, y
      double precision eps, epsj, h, hmax, hmin, pw, save1, save2,
     1                 ss, t, uround, ymax
      double precision d, r, r0, t, yj
      double precision one, rep, zero
      dimension y(n0,1)
c
      common /epcom1/ t,h,hmin,hmax,eps,ss,uround,n,mf,kflag,jstart
      common /epcom2/ ymax(1)
      common /epcom4/ save1(1)
      common /epcom5/ save2(1)
      common /epcom6/ pw(1)
      common /epcom7/ ipiv(1)
      common /epcom8/ epsj,nsq
      data one /1.0d0/, rep /1.0d-3/, zero /0.0d0/
      if (miter .eq. 2) go to 20
c if miter = 1, call pederv and multiply by a scalar.  -----------------
      call pederv (n, t, y, pw, n0)
      do 10 i = 1,nsq
  10    pw(i) = pw(i)*con
      go to 60
c if miter = 2, make n calls to diffun to approximate j. ---------------
  20  d = zero
      do 30 i = 1,n
  30    d = d + save2(i)**2
      r0 = dabs(h)*dsqrt(d)*uround/rep
      j1 = 0
      do 50 j = 1,n
        yj = y(j,1)
        r = epsj*ymax(j)
        r = dmax1(r,r0)
        y(j,1) = y(j,1) + r
        d = con/r
        call diffun (n, t, y, save1)
        do 40 i = 1,n
  40      pw(i+j1) = (save1(i) - save2(i))*d
        y(j,1) = yj
        j1 = j1 + n0
  50    continue
c add on the identity matrix.  -----------------------------------------
  60  j = 1
      do 70 i = 1,n
        pw(j) = pw(j) + one
  70    j = j + (n0 + 1)
c get lu decomposition of p. -------------------------------------------
      call dec (n, n0, pw, ipiv, ier)
      return
c-----------------------  end of subroutine pset -----------------------
      end
      subroutine dec (n, ndim, a, ip, ier)
c-----------------------------------------------------------------------
c matrix triangularization by gaussian elimination.
c input..
c    n = order of matrix.
c    ndim = declared dimension of array  a .
c    a = matrix to be triangularized.
c output..
c    a(i,j), i.le.j = upper triangular factor, u .
c    a(i,j), i.gt.j = multipliers = lower triangular factor, i - l.
c    ip(k), k.lt.n = index of k-th pivot row.
c    ip(n) = (-1)**(number of interchanges) or o .
c    ier = 0 if a nonsingular, or k if a found to be
c          singular at stage k.
c use  sol  to obtain solution of linear system.
c determ(a) = ip(n)*a(1,1)*a(2,2)*...*a(n,n).
c if ip(n)=0, a is singular, sol will divide by zero.
c interchanges finished in u , only partly in l .
c
c reference..
c c. b. moler, algorithm 423, linear equation solver,
c comm. assoc. comput. mach., 15 (1972), p. 274.
c-----------------------------------------------------------------------
      integer ier, ip, n, ndim
      integer i, j, k, kp1, m, nm1
      double precision a
      double precision t
      double precision one, zero
      dimension a(ndim,n),ip(n)
      data one /1.0d0/, zero /0.0d0/
      ier = 0
      ip(n) = 1
      if (n .eq. 1) go to 70
      nm1 = n - 1
      do 60 k = 1,nm1
        kp1 = k + 1
        m = k
        do 10 i = kp1,n
  10      if (dabs(a(i,k)) .gt. dabs(a(m,k))) m = i
        ip(k) = m
        t = a(m,k)
        if (m .eq. k) go to 20
        ip(n) = -ip(n)
        a(m,k) = a(k,k)
        a(k,k) = t
 20     if (t .eq. zero) go to 80
        t = one/t
        do 30 i = kp1,n
 30       a(i,k) = -a(i,k)*t
        do 50 j = kp1,n
          t = a(m,j)
          a(m,j) = a(k,j)
          a(k,j) = t
          if (t .eq. zero) go to 50
          do 40 i = kp1,n
 40         a(i,j) = a(i,j) + a(i,k)*t
 50       continue
 60     continue
 70   k = n
      if (a(n,n) .eq. zero) go to 80
      return
 80   ier = k
      ip(n) = 0
      return
c-----------------------  end of subroutine dec  -----------------------
      end
      subroutine sol (n, ndim, a, b, ip)
c-----------------------------------------------------------------------
c solution of linear system, a*x = b .
c input..
c   n = order of matrix.
c   ndim = declared dimension of array  a .
c   a = triangularized matrix obtained from dec.
c   b = right hand side vector.
c   ip = pivot vector obtained from dec.
c do not use if dec has set ier .ne. 0.
c output..
c   b = solution vector, x .
c-----------------------------------------------------------------------
      integer ip, n, ndim
      integer i, k, kb, km1, kp1, m, nm1
      double precision a, b
      double precision t
      dimension a(ndim, n), b(n), ip(n)
c
      if (n .eq. 1) go to 50
      nm1 = n - 1
      do 20 k = 1,nm1
        kp1 = k + 1
        m = ip(k)
        t = b(m)
        b(m) = b(k)
        b(k) = t
        do 10 i = kp1,n
 10       b(i) = b(i) + a(i,k)*t
 20     continue
      do 40 kb = 1,nm1
        km1 = n - kb
        k = km1 + 1
        b(k) = b(k)/a(k,k)
        t = -b(k)
        do 30 i = 1,km1
 30       b(i) = b(i) + a(i,k)*t
 40     continue
 50   b(1) = b(1)/a(1,1)
      return
c-----------------------  end of subroutine sol  -----------------------
      end