c-----------------------------------------------------------------------
c demonstration program for the lsoibt package.
c this the version of june 22, 1984.
c this version is in double precision.
c
c this program solves a semi-discretized form of the following system
c of three pdes (each similar to a burgers equation)..
c
c   u(i)   =  -(u(1)+u(2)+u(3)) u(i)   +  eta(i) u(i)    (i=1,2,3),
c       t                           x                xx
c
c on the interval  -1 .le. x .le. 1, and with time t .ge. 0.
c the diffusion coefficients are eta(*) = .1, .02, .01.
c the boundary conditions are u(i) = 0 at x = -1 and x = 1 for all i.
c the initial profile for each u(i) is a square wave..
c     u(i) = 0         on 1/2 .lt. abs(x) .le. 1
c     u(i) = amp(i)/2  on abs(x) = 1/2
c     u(i) = amp(i)    on 0 .le. abs(x) .lt. 1/2
c where the amplitudes are amp(*) = .2, .3, .5.
c
c a simplified galerkin treatment of the spatial variable x is used,
c with piecewise linear basis functions on a uniform mesh of 100
c intervals.  the result is a system of odes in the discrete values
c u(i,k) approximating u(i)  (i=1,2,3) at the interior points
c (k = 1,...,99).  the ode-s are..
c
c    .            .        .
c   (u(i,k-1) + 4 u(i,k) + u(i,k+1))/6  =
c
c     -(1/6dx) (c(k-1)dul(i) + 2c(k)(dul(i)+dur(i)) + c(k+1)dur(i))
c
c     + (eta(i)/dx**2) (dur(i) - dul(i))     (i=1,2,3,  k=1,...,99),
c
c where
c     c(j) = u(1,j)+u(2,j)+u(3,j),   dx = .02 = the interval size,
c     dul(i) = u(i,k) - u(i,k-1),   dur(i) = u(i,k+1) - u(i,k).
c terms involving boundary values (subscripts 0 or 100) are dropped
c from the equations for k = 1 and k = 99 above.
c
c the problem is run for each of the 4 values of mf, and for two values
c of the tolerances.  output is taken at t = .1, .2, .3, .4.
c output is on unit lout, set to 6 in a data statement below.
c-----------------------------------------------------------------------
      external res, addabt, jacbt
      integer ncomp, nip, nm1
      integer i, io, istate, itol, iwork, jtol, lout, liw, lrw,
     1   meth, miter, mf, neq, nerr, nint, nout
      double precision eodsq, r6d
      double precision abermx, atol, dx, errfac, eta, hun, one,
     1   rtol, rwork, six, t, tinit, tlast, tout, tols, two, y, ydoti
      dimension eta(3), y(297), ydoti(297), tout(4), tols(2)
      dimension rwork(7447), iwork(317)
c pass problem parameters in the common block par.
      common /par/ r6d, eodsq(3), ncomp, nip, nm1
c
c set problem parameters and run parameters
      data eta/0.1d0,0.02d0,0.01d0/, tinit/0.0d0/, tlast/0.4d0/
      data one/1.0d0/, two/2.0d0/, six/6.0d0/, hun/100.0d0/
      data tout/.10d0,.20d0,.30d0,.40d0/
      data lout/6/, nout/4/, lrw/7447/, liw/317/
      data itol/1/, tols/1.0d-3, 1.0d-6/
c
c set mesh parameters nint, dxc etc.
      nint = 100
      ncomp = 3
      dx = two/dfloat(nint)
      r6d = one/(six*dx)
      do 10 i = 1,ncomp
 10     eodsq(i) = eta(i)/dx**2
      nip = nint - 1
      neq = ncomp*nip
      nm1 = nip - 1
      iwork(1) = ncomp
      iwork(2) = nip
c
      nerr = 0
c
c set the initial conditions (for output purposes only).
      call setic (nint, ncomp, y)
c
      write (lout,1000)
      write (lout,1100) (eta(i),i=1,ncomp), tinit, tlast, nint,
     1   ncomp, nip, neq
      write (lout,1200)
      call edit (y, ncomp, nip, lout)
c
c the jtol loop is over error tolerances.
      do 200 jtol = 1,2
      rtol = tols(jtol)
      atol = rtol
c
c the meth/miter loops cover 4 values of method flag mf.
      do 100 meth = 1,2
       do 100 miter = 1,2
        mf = 10*meth + miter
c
c set the initial conditions.
        call setic (nint, ncomp, y)
        t = tinit
        istate = 0
c
        write (lout,1500)  rtol, atol, mf
c
c loop over output times for each case
        do 80 io = 1,nout
c
          call lsoibt (res, addabt,jacbt, neq, y, ydoti, t, tout(io),
     1     itol,rtol,atol, 1, istate, 0, rwork,lrw,iwork,liw, mf)
c
          write (lout,2000) t, rwork(11), iwork(14), iwork(11)
          if (io .eq. nout) call edit (y, ncomp, nip, lout)
c
c if istate is not 2 on return, go to error section.
          if (istate .ne. 2)  go to 90
c
 80       continue
c
c print final statistics.
        write (lout,3000) mf, iwork(11), iwork(12), iwork(13),
     1                iwork(17), iwork(18)
c
c estimate final error and print result.
        call maxerr (y, ncomp, nip, abermx)
        errfac = abermx/tols(jtol)
        if (errfac .gt. hun)  go to 85
          write (lout,5000) errfac
          go to 100
 85       write (lout,5100) errfac
          nerr = nerr + 1
          go to 100
 90     write (lout,4000) mf, t, istate
        nerr = nerr + 1
 100    continue
 200  continue
c
      write (lout,6000) nerr
      stop
c
 1000 format(1h1/20x,33h demonstration problem for lsoibt//
     1   10x,47hgalerkin method solution of system of 3 pde-s..//
     2   10x,51h  u(i)   =  -(u(1)+u(2)+u(3)) u(i)   +  eta(i) u(i),
     3   5x,9h(i=1,2,3)/16x,1ht,27x,1hx,16x,2hxx//
     4   10x,48hx interval is -1 to 1,  zero boundary conditions/
     5   10x,52hx discretized using piecewise linear basis functions)
 1100 format(/10x,33hfixed parameters are as follows../
     1       13x,32hdiffusion coefficients are eta =,3e10.2/
     2       13x,5ht0 = ,e12.5/13x,8htlast = ,e12.5/
     3       13x,35huniform mesh, number of intervals =,i4/
     4       13x,15hblock size mb =,i2/13x,15hno. blocks nb =,i4/
     5       13x,21hode system size neq =,i5//)
c
 1200 format(/1x,19h initial profiles../)
c
 1500 format(1h1/1x,15hrun with rtol =,e9.1,8h  atol =,e9.1,
     1       7h   mf =,i3///)
c
 2000 format(1x,12hat time t = ,e12.5,16h     current h =,e12.5,
     1       18h   current order =,i2,16h   current nst =,i5/)
c
 3000 format(//1x,26hfinal statistics for mf = ,i2,3h.. ,
     1       i5,7h steps,,i6,5h res,,i6,11h jacobians,,
     2       16h    rwork size =,i6,17h,    iwork size =,i6)
c
 4000 format(//1x,20x,28hfinal time reached for mf = ,i2,
     1       9h was t = ,e12.5/25x,18hat which istate = ,i2//1x,80(1h*))
 5000 format(1x,34hfinal output is correct to within ,e8.1,
     1       30h  times local error tolerance. )
 5100 format(1x,25hfinal output is wrong by ,e8.1,
     1       30h  times local error tolerance. )
 6000 format(////1x,17h run completed.. ,i3,19h errors encountered)
c
c end of main program for the lsoibt demonstration problem.
      end
      subroutine setic (nint, mb, y)
c this routine loads the y array with initial data based on a
c square wave profile for each of the mb pde variables.
c
      integer nint, mb, i, k, nip, n14, n14m1, n14p1, n34, n34m1, n34p1
      double precision y,  amp, half, zero
      dimension y(mb,1), amp(3)
      data zero/0.0d0/, half/0.5d0/, amp/0.2d0,0.3d0,0.5d0/
c
      nip = nint - 1
      n14 = nint/4
      n34 = 3*n14
      n14m1 = n14 - 1
      n14p1 = n14 + 1
      n34m1 = n34 - 1
      n34p1 = n34 + 1
c
      do 15 k = 1,n14m1
        do 10 i = 1,mb
 10       y(i,k) = zero
 15     continue
c
      do 20 i = 1,mb
 20     y(i,n14) = half*amp(i)
c
      do 35 k = n14p1,n34m1
        do 30 i = 1,mb
 30       y(i,k) = amp(i)
 35     continue
c
      do 40 i = 1,mb
 40     y(i,n34) = half*amp(i)
c
      do 55 k = n34p1,nip
        do 50 i = 1,mb
 50       y(i,k) = zero
 55     continue
c
      return
c end of subroutine setic
      end
      subroutine res (n, t, y, v, r, ires)
c this subroutine computes the residual vector
c   r = g(t,y) - a(t,y)*v
c for the lsoibt demonstration problem, using routines gfun and subav.
c if ires = -1, only g(t,y) is returned in r, since a(t,y) does
c not depend on y.
c no changes need to be made to this routine if nip is changed.
c
      integer ires, n,  ncomp, nip, nm1
      double precision t, y, v, r, r6d, eodsq
      dimension y(n), v(n), r(n)
      common /par/ r6d, eodsq(3), ncomp, nip, nm1
c
      call gfun (t, y, r, ncomp)
      if (ires .eq. -1) return
c
      call subav (r, v, ncomp)
c
      return
c end of subroutine res
      end
      subroutine gfun (t, y, g, mb)
c this subroutine computes the function g(y,t) for the
c lsoibt demonstration problem.
c it uses r6d = 1/(6*dx), eodsq(*) = eta(*)/dx**2, nip,
c and nm1 = nip - 1 from the common block par.
c
      integer mb,  ncomp, nip, nm1,  i, k
      double precision t, y, g,  r6d, eodsq,  cc, cl, cr, dli, dri, two
      dimension g(mb,1), y(mb,1)
      common /par/ r6d, eodsq(3), ncomp, nip, nm1
      data two/2.0d0/
c
c left-most interior point (k = 1)
      cc = y(1,1) + y(2,1) + y(3,1)
      cr = y(1,2) + y(2,2) + y(3,2)
      do 10 i = 1,mb
        dri = y(i,2) - y(i,1)
        g(i,1) = -r6d*(two*cc*y(i,2) + cr*dri)
     1         + eodsq(i)*(dri - y(i,1))
 10     continue
c
c interior points k = 2 to nip-1
      do 20 k = 2,nm1
        cl = y(1,k-1) + y(2,k-1) + y(3,k-1)
        cc = y(1,k) + y(2,k) + y(3,k)
        cr = y(1,k+1) + y(2,k+1) + y(3,k+1)
        do 15 i = 1,mb
          dli = y(i,k) - y(i,k-1)
          dri = y(i,k+1) - y(i,k)
          g(i,k) = -r6d*(cl*dli + two*cc*(dli + dri) + cr*dri)
     1           + eodsq(i)*(dri - dli)
 15       continue
 20     continue
c
c right-most interior point (k = nip)
      cl = y(1,nm1) + y(2,nm1) + y(3,nm1)
      cc = y(1,nip) + y(2,nip) + y(3,nip)
      do 30 i = 1,mb
        dli = y(i,nip) - y(i,nm1)
        g(i,nip) = -r6d*(cl*dli - two*cc*y(i,nm1))
     1           - eodsq(i)*(y(i,nip) + dli)
 30     continue
c
      return
c end of subroutine gfun
      end
      subroutine subav (r, v, mb)
c this routine subtracts the matrix a time the vector v from r,
c in order to form the residual vector, stored in r.
c
      integer mb,  ncomp, nip, nm1,  i, k
      double precision r, v,  r6d, eodsq,  aa1, aa4, four, one, six
      dimension r(mb,1), v(mb,1)
      common /par/ r6d, eodsq(3), ncomp, nip, nm1
      data one /1.0d0/, four /4.0d0/, six /6.0d0/
c
      aa1 = one/six
      aa4 = four/six
c
      do 10 i = 1,mb
 10     r(i,1) = r(i,1) - (aa4*v(i,1) + aa1*v(i,2))
c
      do 20 k = 2,nm1
        do 15 i = 1,mb
 15       r(i,k) = r(i,k) - (aa1*v(i,k-1) + aa4*v(i,k) + aa1*v(i,k+1))
 20     continue
c
      do 30 i = 1,mb
 30     r(i,nip) = r(i,nip) - (aa1*v(i,nm1) + aa4*v(i,nip))
c
      return
c end of subroutine subav
      end
      subroutine addabt (n, t, y, mb, nb, pa, pb, pc)
c this subroutine computes the elements of the matrix a,
c and adds them to pa, pb, and pc in the appropriate manner.
c this coding is for the lsoibt demonstration problem.
c the matrix a is tridiagonal, of order n,with
c nonzero elements (reading across) of 1/6,4/6,1/6.
c
      integer n, mb, nb,  i, k
      double precision pa, pb, pc, t, y,  aa1, aa4, four, one, six
      dimension y(mb,nb),pa(mb,mb,nb),pb(mb,mb,nb),pc(mb,mb,nb)
      data one/1.0d0/, four/4.0d0/, six/6.0d0/
c
      aa1 = one/six
      aa4 = four/six
      do 50 k = 1,nb
        do 10 i = 1,mb
 10       pa(i,i,k) = pa(i,i,k) + aa4
        if (k .eq. nb) go to 25
        do 20 i = 1,mb
 20       pb(i,i,k) = pb(i,i,k) + aa1
 25     continue
        if (k .eq. 1) go to 50
        do 30 i = 1,mb
 30       pc(i,i,k) = pc(i,i,k) + aa1
 50     continue
c
      return
c end of subroutine addabt
      end
      subroutine jacbt (n, t, y, s, mb, nb, pa, pb, pc)
c this subroutine computes the jacobian dg/dy=d(g-a*s)/dy
c which has block tridiagonal structure. the main, upper,
c and lower diagonals are stored in pa, pb, and pc respectively.
c
      integer n, mb, nb,  ncomp, nip, nm1,  i, j, k
      double precision t, y, s, pa, pb, pc,  r6d, eodsq,  cc, cl, cr,
     1   dlj, drj, paij, pbij, pcij, terma, termb, termc, two
      dimension y(mb,nb),s(n),pa(mb,mb,nb),pb(mb,mb,nb),pc(mb,mb,nb)
      common /par/ r6d, eodsq(3), ncomp, nip, nm1
      data two/2.0d0/
c
c left-most interior point (k = 1)
      cc = y(1,1) + y(2,1) + y(3,1)
      cr = y(1,2) + y(2,2) + y(3,2)
      terma = r6d*cr
      termb = -r6d*(two*cc + cr)
      do 20 j = 1,mb
        drj = y(j,2) - y(j,1)
        paij = -r6d*two*y(j,2)
        pbij = -r6d*drj
        do 10 i = 1,mb
          pa(i,j,1) = paij
 10       pb(i,j,1) = pbij
        pa(j,j,1) = pa(j,j,1) + terma - two*eodsq(j)
        pb(j,j,1) = pb(j,j,1) + termb + eodsq(j)
 20     continue
c
c interior points k = 2 to nip-1
      do 50 k = 2,nm1
        cl = y(1,k-1) + y(2,k-1) + y(3,k-1)
        cc = y(1,k) + y(2,k) + y(3,k)
        cr = y(1,k+1) + y(2,k+1) + y(3,k+1)
        terma = r6d*(cr - cl)
        termb = -r6d*(two*cc + cr)
        termc = r6d*(two*cc + cl)
        do 40 j = 1,mb
          dlj = y(j,k) - y(j,k-1)
          drj = y(j,k+1) - y(j,k)
          paij = -r6d*two*(dlj + drj)
          pbij = -r6d*drj
          pcij = -r6d*dlj
          do 30 i = 1,mb
            pa(i,j,k) = paij
            pb(i,j,k) = pbij
 30         pc(i,j,k) = pcij
          pa(j,j,k) = pa(j,j,k) + terma - two*eodsq(j)
          pb(j,j,k) = pb(j,j,k) + termb + eodsq(j)
          pc(j,j,k) = pc(j,j,k) + termc + eodsq(j)
 40       continue
 50     continue
c
c right-most interior point (k = nip)
      cl = y(1,nm1) + y(2,nm1) + y(3,nm1)
      cc = y(1,nip) + y(2,nip) + y(3,nip)
      terma = -r6d*cl
      termc = r6d*(two*cc + cl)
      do 70 j = 1,mb
        dlj = y(j,nip) - y(j,nm1)
        paij = r6d*two*y(j,nm1)
        pcij = -r6d*dlj
        do 60 i = 1,mb
          pa(i,j,nip) = paij
 60       pc(i,j,nip) = pcij
        pa(j,j,nip) = pa(j,j,nip) + terma - two*eodsq(j)
        pc(j,j,nip) = pc(j,j,nip) + termc + eodsq(j)
 70     continue
c
      return
c end of subroutine jacbt
      end
      subroutine edit (y, mb, nip, lout)
c this routine prints output.  for each of the mb pde components,
c the values at the nip points are printed.
c all output is on unit lout.
c
      integer mb, nip, lout,  i, k
      double precision y
      dimension y(mb,nip)
c
      do 10 i = 1,mb
        write (lout,20) i, (y(i,k),k=1,nip)
 10     continue
c
 20   format(1x,28h values of pde component i =,i3/10(10e12.4/))
c
      return
c end of subroutine edit
      end
      subroutine maxerr (y, mb, nb, abermx)
c this routine computes the maximum absolute error in the
c y array, as a computed solution at t = 0.4, using data-loaded values
c for accurate answers (from running with much smaller tolerances).
c
      integer mb, nb,  k
      double precision y, abermx,  ae1, ae2, ae3, u1, u2, u3, zero,
     1   u1a, u1b, u1c, u1d, u1e, u1f, u1g,
     2   u2a, u2b, u2c, u2d, u2e, u2f, u2g,
     3   u3a, u3b, u3c, u3d, u3e, u3f, u3g
      dimension y(mb,nb), u1(99), u2(99), u3(99)
      dimension u1a(16),u1b(16),u1c(16),u1d(16),u1e(16),u1f(16),u1g(3),
     2          u2a(16),u2b(16),u2c(16),u2d(16),u2e(16),u2f(16),u2g(3),
     3          u3a(16),u3b(16),u3c(16),u3d(16),u3e(16),u3f(16),u3g(3)
      equivalence (u1a(1),u1(1)), (u1b(1),u1(17)),
     1      (u1c(1),u1(33)), (u1d(1),u1(49)), (u1e(1),u1(65)),
     1      (u1f(1),u1(81)), (u1g(1),u1(97))
      equivalence (u2a(1),u2(1)), (u2b(1),u2(17)),
     1      (u2c(1),u2(33)), (u2d(1),u2(49)), (u2e(1),u2(65)),
     1      (u2f(1),u2(81)), (u2g(1),u2(97))
      equivalence (u3a(1),u3(1)), (u3b(1),u3(17)),
     1      (u3c(1),u3(33)), (u3d(1),u3(49)), (u3e(1),u3(65)),
     1      (u3f(1),u3(81)), (u3g(1),u3(97))
c
      data u1a /
     1  1.70956682d-03, 3.43398445d-03, 5.18783349d-03, 6.98515842d-03,
     1  8.83921016d-03, 1.07622016d-02, 1.27650806d-02, 1.48573251d-02,
     1  1.70467655d-02, 1.93394396d-02, 2.17394852d-02, 2.42490773d-02,
     1  2.68684152d-02, 2.95957660d-02, 3.24275691d-02, 3.53586054d-02/
      data u1b /
     1  3.83822285d-02, 4.14906520d-02, 4.46752791d-02, 4.79270545d-02,
     1  5.12368132d-02, 5.45956048d-02, 5.79949684d-02, 6.14271460d-02,
     1  6.48852271d-02, 6.83632267d-02, 7.18561029d-02, 7.53597274d-02,
     1  7.88708192d-02, 8.23868545d-02, 8.59059616d-02, 8.94268082d-02/
      data u1c /
     1  9.29484864d-02, 9.64703968d-02, 9.99921344d-02, 1.03513375d-01,
     1  1.07033760d-01, 1.10552783d-01, 1.14069668d-01, 1.17583246d-01,
     1  1.21091827d-01, 1.24593066d-01, 1.28083828d-01, 1.31560049d-01,
     1  1.35016617d-01, 1.38447256d-01, 1.41844451d-01, 1.45199401d-01/
      data u1d /
     1  1.48502033d-01, 1.51741065d-01, 1.54904135d-01, 1.57977973d-01,
     1  1.60948623d-01, 1.63801670d-01, 1.66522463d-01, 1.69096305d-01,
     1  1.71508595d-01, 1.73744902d-01, 1.75790974d-01, 1.77632682d-01,
     1  1.79255895d-01, 1.80646319d-01, 1.81789276d-01, 1.82669470d-01/
      data u1e /
     1  1.83270725d-01, 1.83575716d-01, 1.83565712d-01, 1.83220322d-01,
     1  1.82517279d-01, 1.81432251d-01, 1.79938706d-01, 1.78007835d-01,
     1  1.75608540d-01, 1.72707519d-01, 1.69269456d-01, 1.65257378d-01,
     1  1.60633244d-01, 1.55358941d-01, 1.49398029d-01, 1.42718981d-01/
      data u1f /
     1  1.35301474d-01, 1.27148627d-01, 1.18308730d-01, 1.08905085d-01,
     1  9.91559295d-02, 8.93515884d-02, 7.97824293d-02, 7.06663514d-02,
     1  6.21244732d-02, 5.41994827d-02, 4.68848207d-02, 4.01465202d-02,
     1  3.39357642d-02, 2.81954415d-02, 2.28635569d-02, 1.78750916d-02/
      data u1g / 1.31630892d-02, 8.65933391d-03, 4.29480447d-03/
      data u2a /
     1  7.17416019d-06, 1.70782645d-05, 3.31245126d-05, 6.01588363d-05,
     1  1.05339286d-04, 1.79174771d-04, 2.96719122d-04, 4.78862606d-04,
     1  7.53598916d-04, 1.15707860d-03, 1.73420412d-03, 2.53849668d-03,
     1  3.63099110d-03, 5.07800919d-03, 6.94782549d-03, 9.30645443d-03/
      data u2b /
     1  1.22130079d-02, 1.57152366d-02, 1.98459102d-02, 2.46205841d-02,
     1  3.00370492d-02, 3.60764461d-02, 4.27057301d-02, 4.98809820d-02,
     1  5.75510102d-02, 6.56607602d-02, 7.41541974d-02, 8.29764928d-02,
     1  9.20754824d-02, 1.01402468d-01, 1.10912474d-01, 1.20564094d-01/
      data u2c /
     1  1.30319039d-01, 1.40141489d-01, 1.49997326d-01, 1.59853293d-01,
     1  1.69676126d-01, 1.79431680d-01, 1.89084097d-01, 1.98595037d-01,
     1  2.07923034d-01, 2.17023055d-01, 2.25846345d-01, 2.34340694d-01,
     1  2.42451240d-01, 2.50121934d-01, 2.57297724d-01, 2.63927433d-01/
      data u2d /
     1  2.69967170d-01, 2.75383917d-01, 2.80158840d-01, 2.84289739d-01,
     1  2.87792167d-01, 2.90698875d-01, 2.93057586d-01, 2.94927384d-01,
     1  2.96374262d-01, 2.97466488d-01, 2.98270390d-01, 2.98847025d-01,
     1  2.99249945d-01, 2.99524080d-01, 2.99705593d-01, 2.99822450d-01/
      data u2e /
     1  2.99895431d-01, 2.99939301d-01, 2.99963931d-01, 2.99975129d-01,
     1  2.99974996d-01, 2.99961526d-01, 2.99927041d-01, 2.99854809d-01,
     1  2.99712769d-01, 2.99442742d-01, 2.98942676d-01, 2.98038511d-01,
     1  2.96441259d-01, 2.93684573d-01, 2.89040478d-01, 2.81421884d-01/
      data u2f /
     1  2.69315148d-01, 2.50874185d-01, 2.24457680d-01, 1.89885662d-01,
     1  1.49894358d-01, 1.09927672d-01, 7.54041273d-02, 4.90259517d-02,
     1  3.06080023d-02, 1.85165524d-02, 1.09104125d-02, 6.27726960d-03,
     1  3.53002680d-03, 1.94049735d-03, 1.04218859d-03, 5.45964314d-04/
      data u2g / 2.77379128d-04, 1.33343739d-04, 5.32660444d-05/
      data u3a /
     1  1.86765383d-10, 1.96772458d-09, 1.19111389d-08, 5.54964761d-08,
     1  2.18340713d-07, 7.55899524d-07, 2.35604385d-06, 6.70801745d-06,
     1  1.76224112d-05, 4.30351929d-05, 9.82592148d-05, 2.10736217d-04,
     1  4.26209304d-04, 8.15657041d-04, 1.48160943d-03, 2.56186555d-03/
      data u3b /
     1  4.22851247d-03, 6.68078970d-03, 1.01317466d-02, 1.47903961d-02,
     1  2.08424987d-02, 2.84336008d-02, 3.76573037d-02, 4.85502549d-02,
     1  6.10936693d-02, 7.52198901d-02, 9.08218891d-02, 1.07763660d-01,
     1  1.25889931d-01, 1.45034247d-01, 1.65025016d-01, 1.85689556d-01/
      data u3c /
     1  2.06856371d-01, 2.28356037d-01, 2.50021072d-01, 2.71685149d-01,
     1  2.93181998d-01, 3.14344301d-01, 3.35002907d-01, 3.54986687d-01,
     1  3.74123404d-01, 3.92241969d-01, 4.09176451d-01, 4.24772089d-01,
     1  4.38893320d-01, 4.51433444d-01, 4.62324969d-01, 4.71549073d-01/
      data u3d /
     1  4.79142163d-01, 4.85197409d-01, 4.89859810d-01, 4.93314543d-01,
     1  4.95770115d-01, 4.97439231d-01, 4.98520996d-01, 4.99187563d-01,
     1  4.99576941d-01, 4.99791928d-01, 4.99903753d-01, 4.99958343d-01,
     1  4.99983239d-01, 4.99993785d-01, 4.99997902d-01, 4.99999367d-01/
      data u3e /
     1  4.99999835d-01, 4.99999965d-01, 4.99999995d-01, 5.00000000d-01,
     1  5.00000000d-01, 4.99999997d-01, 4.99999976d-01, 4.99999863d-01,
     1  4.99999315d-01, 4.99996914d-01, 4.99987300d-01, 4.99951740d-01,
     1  4.99829328d-01, 4.99435130d-01, 4.98245007d-01, 4.94883400d-01/
      data u3f /
     1  4.86081966d-01, 4.65174923d-01, 4.21856650d-01, 3.47885738d-01,
     1  2.49649938d-01, 1.51648615d-01, 7.80173239d-02, 3.47983164d-02,
     1  1.38686441d-02, 5.05765688d-03, 1.71052539d-03, 5.38966324d-04,
     1  1.57923694d-04, 4.27352191d-05, 1.05512005d-05, 2.33068621d-06/
      data u3g / 4.45404604d-07, 6.88336884d-08, 7.23875975d-09/
      data zero/0.0d0/
c
      abermx = zero
      do 10 k = 1,99
        ae1 = dabs(y(1,k) - u1(k))
        ae2 = dabs(y(2,k) - u2(k))
        ae3 = dabs(y(3,k) - u3(k))
        abermx = dmax1(abermx, ae1, ae2, ae3)
 10     continue
c
c
      return
c end of subroutine maxerr
c end of demonstration program for lsoibt.
      end