share/odepack/fortran/opkdemo8.f

   1       program opkdemo8
   2 c-----------------------------------------------------------------------
   3 c Demonstration program for the DLSOIBT package.
   4 c This is the version of 14 June 2001.
   5 c
   6 c This version is in double precision.
   7 c
   8 c This program solves a semi-discretized form of the following system
   9 c Of three PDEs (each similar to a Burgers equation):
  10 c
  11 c   u(i)   =  -(u(1)+u(2)+u(3)) u(i)   +  eta(i) u(i)    (i=1,2,3),
  12 c       t                           x                xx
  13 c
  14 c on the interval  -1 .le. x .le. 1, and with time t .ge. 0.
  15 c The diffusion coefficients are eta(*) = .1, .02, .01.
  16 c The boundary conditions are u(i) = 0 at x = -1 and x = 1 for all i.
  17 c The initial profile for each u(i) is a square wave:
  18 c     u(i) = 0         on 1/2 .lt. abs(x) .le. 1
  19 c     u(i) = amp(i)/2  on abs(x) = 1/2
  20 c     u(i) = amp(i)    on 0 .le. abs(x) .lt. 1/2
  21 c where the amplitudes are amp(*) = .2, .3, .5.
  22 c
  23 c A simplified Galerkin treatment of the spatial variable x is used,
  24 c with piecewise linear basis functions on a uniform mesh of 100
  25 c intervals.  The result is a system of ODEs in the discrete values
  26 c u(i,k) approximating u(i)  (i=1,2,3) at the interior points
  27 c (k = 1,...,99).  The ODEs are:
  28 c
  29 c    .            .        .
  30 c   (u(i,k-1) + 4 u(i,k) + u(i,k+1))/6  =
  31 c
  32 c     -(1/6dx) (c(k-1)dul(i) + 2c(k)(dul(i)+dur(i)) + c(k+1)dur(i))
  33 c
  34 c     + (eta(i)/dx**2) (dur(i) - dul(i))     (i=1,2,3,  k=1,...,99),
  35 c
  36 c where
  37 c     c(j) = u(1,j)+u(2,j)+u(3,j),   dx = .02 = the interval size,
  38 c     dul(i) = u(i,k) - u(i,k-1),   dur(i) = u(i,k+1) - u(i,k).
  39 c Terms involving boundary values (subscripts 0 or 100) are dropped
  40 c from the equations for k = 1 and k = 99 above.
  41 c
  42 c The problem is run for each of the 4 values of mf, and for two values
  43 c of the tolerances.  Output is taken at t = .1, .2, .3, .4.
  44 c Output is on unit lout, set to 6 in a data statement below.
  45 c-----------------------------------------------------------------------
  46       external res, addabt, jacbt
  47       integer ncomp, nip, nm1
  48       integer i, io, istate, itol, iwork, jtol, lout, liw, lrw,
  49      1   meth, miter, mf, neq, nerr, nint, nout
  50       double precision eodsq, r6d
  51       double precision abermx, atol, dx, errfac, eta, hun, one,
  52      1   rtol, rwork, six, t, tinit, tlast, tout, tols, two, y, ydoti
  53       dimension eta(3), y(297), ydoti(297), tout(4), tols(2)
  54       dimension rwork(7447), iwork(317)
  55       dimension neq(1), rtol(1), atol(1)
  56 c Pass problem parameters in the common block par.
  57       common /par/ r6d, eodsq(3), ncomp, nip, nm1
  58 c
  59 c Set problem parameters and run parameters
  60       data eta/0.1d0,0.02d0,0.01d0/, tinit/0.0d0/, tlast/0.4d0/
  61       data one/1.0d0/, two/2.0d0/, six/6.0d0/, hun/100.0d0/
  62       data tout/.10d0,.20d0,.30d0,.40d0/
  63       data lout/6/, nout/4/, lrw/7447/, liw/317/
  64       data itol/1/, tols/1.0d-3, 1.0d-6/
  65 c
  66 c Set mesh parameters nint, dxc etc.
  67       nint = 100
  68       ncomp = 3
  69       dx = two/nint
  70       r6d = one/(six*dx)
  71       do 10 i = 1,ncomp
  72  10     eodsq(i) = eta(i)/dx**2
  73       nip = nint - 1
  74       neq(1) = ncomp*nip
  75       nm1 = nip - 1
  76       iwork(1) = ncomp
  77       iwork(2) = nip
  78 c
  79       nerr = 0
  80 c
  81 c Set the initial conditions (for output purposes only).
  82       call setic (nint, ncomp, y)
  83 c
  84       write (lout,1000)
  85       write (lout,1100) (eta(i),i=1,ncomp), tinit, tlast, nint,
  86      1      ncomp, nip, neq(1)
  87       write (lout,1200)
  88       call edit (y, ncomp, nip, lout)
  89 c
  90 c The jtol loop is over error tolerances.
  91       do 200 jtol = 1,2
  92       rtol(1) = tols(jtol)
  93       atol(1) = rtol(1)
  94 c
  95 c The meth/miter loops cover 4 values of method flag mf.
  96       do 100 meth = 1,2
  97        do 100 miter = 1,2
  98         mf = 10*meth + miter
  99 c
 100 c Set the initial conditions.
 101         call setic (nint, ncomp, y)
 102         t = tinit
 103         istate = 0
 104 c
 105         write (lout,1500)  rtol(1), atol(1), mf
 106 c
 107 c Loop over output times for each case
 108         do 80 io = 1,nout
 109 c
 110           call dlsoibt (res, addabt,jacbt, neq, y, ydoti, t, tout(io),
 111      1     itol,rtol,atol, 1, istate, 0, rwork,lrw,iwork,liw, mf)
 112 c
 113           write (lout,2000) t, rwork(11), iwork(14), iwork(11)
 114           if (io .eq. nout) call edit (y, ncomp, nip, lout)
 115 c
 116 c If istate is not 2 on return, print message and go to next case.
 117           if (istate .ne. 2) then
 118             write (lout,4000) mf, t, istate
 119             nerr = nerr + 1
 120             go to 100
 121             endif
 122 c
 123  80       continue
 124 c
 125 c Print final statistics.
 126         write (lout,3000) mf, iwork(11), iwork(12), iwork(13),
 127      1         iwork(17), iwork(18)
 128 c
 129 c Estimate final error and print result.
 130         call maxerr (y, ncomp, nip, abermx)
 131         errfac = abermx/tols(jtol)
 132         if (errfac .lt. hun) then
 133           write (lout,5000) errfac
 134         else
 135           write (lout,5100) errfac
 136           nerr = nerr + 1
 137           endif
 138  100    continue
 139  200  continue
 140 c
 141       write (lout,6000) nerr
 142 c      stop
 143 c
 144  1000 format(/20x,' Demonstration Problem for DLSOIBT'//
 145      1   10x,'Galerkin method solution of system of 3 PDEs:'//
 146      2   10x,'  u(i)   =  -(u(1)+u(2)+u(3)) u(i)   +  eta(i) u(i)',
 147      3   5x,'(i=1,2,3)',/16x,'t',27x,'x',16x,'xx'//
 148      4   10x,'x interval is -1 to 1,  zero boundary conditions'/
 149      5   10x,'x discretized using piecewise linear basis functions')
 150  1100 format(/10x,'Fixed parameters are as follows:'/
 151      1       13x,'Diffusion coefficients are eta =',3d10.2/
 152      2       13x,'t0 = ',d12.5/13x,'tlast = ',d12.5/
 153      3       13x,'Uniform mesh, number of intervals =',i4/
 154      4       13x,'Block size mb =',i2/13x,'Number of blocks nb =',i4/
 155      5       13x,'ODE system size neq =',i5//)
 156 c
 157  1200 format(/'Initial profiles:'/)
 158 c
 159  1500 format(////90('*')//'Run with rtol =',d9.1,'  atol =',d9.1,
 160      1       '   mf =',i3///)
 161 c
 162  2000 format(' At time t =',d12.5,'  current h =',d12.5,
 163      1       '  current order =',i2,'  current nst =',i5/)
 164 c
 165  3000 format(//'Final statistics for mf = ',i2,':',
 166      1       i5,' steps,',i6,' res,',i6,' jacobians,'/
 167      2       30x,       'rwork size =',i6,',  iwork size =',i6)
 168 c
 169  4000 format(//20x,'Final time reached for mf = ',i2,
 170      1       ' was t = ',d12.5/25x,'at which istate = ',i2//)
 171  5000 format('Final output is correct to within ',d9.2,
 172      1       '  times local error tolerance. ')
 173  5100 format('Final output is wrong by ',d9.2,
 174      1       '  times local error tolerance.')
 175  6000 format(//90('*')//'Run completed: ',i3,' errors encountered')
 176 c
 177 c end of main program for the DLSOIBT demonstration problem
 178       end
 179
 180       subroutine setic (nint, mb, y)
 181 c This routine loads the y array with initial data based on a
 182 c square wave profile for each of the mb PDE variables.
 183 c
 184       integer nint, mb, i, k, nip, n14, n34
 185       double precision y,  amp, half, zero
 186       dimension y(mb,*), amp(3)
 187       data zero/0.0d0/, half/0.5d0/, amp/0.2d0,0.3d0,0.5d0/
 188 c
 189       nip = nint - 1
 190       n14 = nint/4
 191       n34 = 3*n14
 192 c
 193       do 15 k = 1,n14-1
 194         do 10 i = 1,mb
 195  10       y(i,k) = zero
 196  15     continue
 197 c
 198       do 20 i = 1,mb
 199  20     y(i,n14) = half*amp(i)
 200 c
 201       do 35 k = n14+1,n34-1
 202         do 30 i = 1,mb
 203  30       y(i,k) = amp(i)
 204  35     continue
 205 c
 206       do 40 i = 1,mb
 207  40     y(i,n34) = half*amp(i)
 208 c
 209       do 55 k = n34+1,nip
 210         do 50 i = 1,mb
 211  50       y(i,k) = zero
 212  55     continue
 213 c
 214       return
 215 c end of subroutine setic
 216       end
 217
 218       subroutine res (n, t, y, v, r, ires)
 219 c This subroutine computes the residual vector
 220 c   r = g(t,y) - A(t,y)*v
 221 c using routines gfun and subav.
 222 c If ires = -1, only g(t,y) is returned in r, since A(t,y) does
 223 c not depend on y.
 224 c No changes need to be made to this routine if nip is changed.
 225 c
 226       integer ires, n,  ncomp, nip, nm1
 227       double precision t, y, v, r, r6d, eodsq
 228       dimension y(*), v(*), r(*)
 229       dimension n(*)
 230       common /par/ r6d, eodsq(3), ncomp, nip, nm1
 231 c
 232       call gfun (t, y, r, ncomp)
 233       if (ires .eq. -1) return
 234 c
 235       call subav (r, v, ncomp)
 236 c
 237       return
 238 c end of subroutine res
 239       end
 240
 241       subroutine gfun (t, y, g, mb)
 242 c This subroutine computes the right-hand side function g(y,t).
 243 c It uses r6d = 1/(6*dx), eodsq(*) = eta(*)/dx**2, nip,
 244 c and nm1 = nip - 1 from the Common block par.
 245 c
 246       integer mb,  ncomp, nip, nm1,  i, k
 247       double precision t, y, g,  r6d, eodsq,  cc, cl, cr, dli, dri, two
 248       dimension g(mb,*), y(mb,*)
 249       common /par/ r6d, eodsq(3), ncomp, nip, nm1
 250       data two/2.0d0/
 251 c
 252 c left-most interior point (k = 1)
 253       cc = y(1,1) + y(2,1) + y(3,1)
 254       cr = y(1,2) + y(2,2) + y(3,2)
 255       do 10 i = 1,mb
 256         dri = y(i,2) - y(i,1)
 257         g(i,1) = -r6d*(two*cc*y(i,2) + cr*dri)
 258      1         + eodsq(i)*(dri - y(i,1))
 259  10     continue
 260 c
 261 c interior points k = 2 to nip-1
 262       do 20 k = 2,nm1
 263         cl = y(1,k-1) + y(2,k-1) + y(3,k-1)
 264         cc = y(1,k) + y(2,k) + y(3,k)
 265         cr = y(1,k+1) + y(2,k+1) + y(3,k+1)
 266         do 15 i = 1,mb
 267           dli = y(i,k) - y(i,k-1)
 268           dri = y(i,k+1) - y(i,k)
 269           g(i,k) = -r6d*(cl*dli + two*cc*(dli + dri) + cr*dri)
 270      1           + eodsq(i)*(dri - dli)
 271  15       continue
 272  20     continue
 273 c
 274 c right-most interior point (k = nip)
 275       cl = y(1,nm1) + y(2,nm1) + y(3,nm1)
 276       cc = y(1,nip) + y(2,nip) + y(3,nip)
 277       do 30 i = 1,mb
 278         dli = y(i,nip) - y(i,nm1)
 279         g(i,nip) = -r6d*(cl*dli - two*cc*y(i,nm1))
 280      1           - eodsq(i)*(y(i,nip) + dli)
 281  30     continue
 282 c
 283         t = t
 284       return
 285 c end of subroutine gfun
 286       end
 287
 288       subroutine subav (r, v, mb)
 289 c This routine subtracts the matrix a time the vector v from r,
 290 c in order to form the residual vector, stored in r.
 291 c
 292       integer mb,  ncomp, nip, nm1,  i, k
 293       double precision r, v,  r6d, eodsq,  aa1, aa4, four, one, six
 294       dimension r(mb,*), v(mb,*)
 295       common /par/ r6d, eodsq(3), ncomp, nip, nm1
 296       data one /1.0d0/, four /4.0d0/, six /6.0d0/
 297 c
 298       aa1 = one/six
 299       aa4 = four/six
 300 c
 301       do 10 i = 1,mb
 302  10     r(i,1) = r(i,1) - (aa4*v(i,1) + aa1*v(i,2))
 303 c
 304       do 20 k = 2,nm1
 305         do 15 i = 1,mb
 306  15       r(i,k) = r(i,k) - (aa1*v(i,k-1) + aa4*v(i,k) + aa1*v(i,k+1))
 307  20     continue
 308 c
 309       do 30 i = 1,mb
 310  30     r(i,nip) = r(i,nip) - (aa1*v(i,nm1) + aa4*v(i,nip))
 311 c
 312       return
 313 c end of subroutine subav
 314       end
 315
 316       subroutine addabt (n, t, y, mb, nb, pa, pb, pc)
 317 c This subroutine computes the elements of the matrix A,
 318 c and adds them to pa, pb, and pc in the appropriate manner.
 319 c The matrix A is tridiagonal, of order n, with
 320 c nonzero elements (reading across) of 1/6, 4/6, 1/6.
 321 c
 322       integer n, mb, nb, i, k
 323       double precision pa, pb, pc, t, y,  aa1, aa4, four, one, six
 324       dimension y(mb,nb),pa(mb,mb,nb),pb(mb,mb,nb),pc(mb,mb,nb)
 325       data one/1.0d0/, four/4.0d0/, six/6.0d0/
 326 c
 327       aa1 = one/six
 328       aa4 = four/six
 329       do 50 k = 1,nb
 330         do 10 i = 1,mb
 331  10       pa(i,i,k) = pa(i,i,k) + aa4
 332         if (k .ne. nb) then
 333           do 20 i = 1,mb
 334  20         pb(i,i,k) = pb(i,i,k) + aa1
 335         endif
 336         if (k .ne. 1) then
 337           do 30 i = 1,mb
 338  30         pc(i,i,k) = pc(i,i,k) + aa1
 339         endif
 340  50     continue
 341 c
 342       return
 343 c end of subroutine addabt
 344       end
 345
 346       subroutine jacbt (n, t, y, s, mb, nb, pa, pb, pc)
 347 c This subroutine computes the Jacobian dg/dy = d(g-A*s)/dy
 348 c which has block-tridiagonal structure.  The main, upper, and
 349 c lower diagonals are stored in pa, pb, and pc, respectively.
 350 c
 351       integer n, mb, nb,  ncomp, nip, nm1, i, j, k
 352       double precision t, y, s, pa, pb, pc,  r6d, eodsq,  cc, cl, cr,
 353      1   dlj, drj, paij, pbij, pcij, terma, termb, termc, two
 354       dimension y(mb,nb),s(n),pa(mb,mb,nb),pb(mb,mb,nb),pc(mb,mb,nb)
 355       common /par/ r6d, eodsq(3), ncomp, nip, nm1
 356       data two/2.0d0/
 357 c
 358 c left-most interior point (k = 1)
 359       cc = y(1,1) + y(2,1) + y(3,1)
 360       cr = y(1,2) + y(2,2) + y(3,2)
 361       terma = r6d*cr
 362       termb = -r6d*(two*cc + cr)
 363       do 20 j = 1,mb
 364         drj = y(j,2) - y(j,1)
 365         paij = -r6d*two*y(j,2)
 366         pbij = -r6d*drj
 367         do 10 i = 1,mb
 368           pa(i,j,1) = paij
 369  10       pb(i,j,1) = pbij
 370         pa(j,j,1) = pa(j,j,1) + terma - two*eodsq(j)
 371         pb(j,j,1) = pb(j,j,1) + termb + eodsq(j)
 372  20     continue
 373 c
 374 c interior points k = 2 to nip-1
 375       do 50 k = 2,nm1
 376         cl = y(1,k-1) + y(2,k-1) + y(3,k-1)
 377         cc = y(1,k) + y(2,k) + y(3,k)
 378         cr = y(1,k+1) + y(2,k+1) + y(3,k+1)
 379         terma = r6d*(cr - cl)
 380         termb = -r6d*(two*cc + cr)
 381         termc = r6d*(two*cc + cl)
 382         do 40 j = 1,mb
 383           dlj = y(j,k) - y(j,k-1)
 384           drj = y(j,k+1) - y(j,k)
 385           paij = -r6d*two*(dlj + drj)
 386           pbij = -r6d*drj
 387           pcij = -r6d*dlj
 388           do 30 i = 1,mb
 389             pa(i,j,k) = paij
 390             pb(i,j,k) = pbij
 391  30         pc(i,j,k) = pcij
 392           pa(j,j,k) = pa(j,j,k) + terma - two*eodsq(j)
 393           pb(j,j,k) = pb(j,j,k) + termb + eodsq(j)
 394           pc(j,j,k) = pc(j,j,k) + termc + eodsq(j)
 395  40       continue
 396  50     continue
 397 c
 398 c right-most interior point (k = nip)
 399       cl = y(1,nm1) + y(2,nm1) + y(3,nm1)
 400       cc = y(1,nip) + y(2,nip) + y(3,nip)
 401       terma = -r6d*cl
 402       termc = r6d*(two*cc + cl)
 403       do 70 j = 1,mb
 404         dlj = y(j,nip) - y(j,nm1)
 405         paij = r6d*two*y(j,nm1)
 406         pcij = -r6d*dlj
 407         do 60 i = 1,mb
 408           pa(i,j,nip) = paij
 409  60       pc(i,j,nip) = pcij
 410         pa(j,j,nip) = pa(j,j,nip) + terma - two*eodsq(j)
 411         pc(j,j,nip) = pc(j,j,nip) + termc + eodsq(j)
 412  70     continue
 413 c
 414       return
 415 c end of subroutine jacbt
 416       end
 417
 418       subroutine edit (y, mb, nip, lout)
 419 c This routine prints output.  For each of the mb PDE components, the
 420 c values at the nip points are printed.  All output is on unit lout.
 421 c
 422       integer mb, nip, lout,  i, k
 423       double precision y
 424       dimension y(mb,nip)
 425 c
 426       do 10 i = 1,mb
 427  10      write (lout,20) i, (y(i,k),k=1,nip)
 428 c
 429  20   format(' Values of PDE component i =',i3/15(7d12.4/))
 430 c
 431       return
 432 c end of subroutine edit
 433       end
 434
 435       subroutine maxerr (y, mb, nb, abermx)
 436 c This routine computes the maximum absolute error in the y array,
 437 c as a computed solution at t = 0.4, using data-loaded values for
 438 c accurate answers (from running with much smaller tolerances).
 439 c
 440       integer mb, nb,  k
 441       double precision y, abermx,  ae1, ae2, ae3, u1, u2, u3, zero
 442 c     1   u1a, u1b, u1c, u1d, u1e, u1f, u1g,
 443 c     2   u2a, u2b, u2c, u2d, u2e, u2f, u2g,
 444 c     3   u3a, u3b, u3c, u3d, u3e, u3f, u3g
 445       dimension y(mb,nb), u1(99), u2(99), u3(99)
 446 c      dimension u1a(16),u1b(16),u1c(16),u1d(16),u1e(16),u1f(16),u1g(3),
 447 c     2          u2a(16),u2b(16),u2c(16),u2d(16),u2e(16),u2f(16),u2g(3),
 448 c     3          u3a(16),u3b(16),u3c(16),u3d(16),u3e(16),u3f(16),u3g(3)
 449 c      equivalence (u1a(1),u1(1)), (u1b(1),u1(17)),
 450 c     1      (u1c(1),u1(33)), (u1d(1),u1(49)), (u1e(1),u1(65)),
 451 c     1      (u1f(1),u1(81)), (u1g(1),u1(97))
 452 c      equivalence (u2a(1),u2(1)), (u2b(1),u2(17)),
 453 c     1      (u2c(1),u2(33)), (u2d(1),u2(49)), (u2e(1),u2(65)),
 454 c     1      (u2f(1),u2(81)), (u2g(1),u2(97))
 455 c      equivalence (u3a(1),u3(1)), (u3b(1),u3(17)),
 456 c     1      (u3c(1),u3(33)), (u3d(1),u3(49)), (u3e(1),u3(65)),
 457 c     1      (u3f(1),u3(81)), (u3g(1),u3(97))
 458 c
 459       data u1 /
 460      1  1.70956682d-03, 3.43398445d-03, 5.18783349d-03, 6.98515842d-03,
 461      1  8.83921016d-03, 1.07622016d-02, 1.27650806d-02, 1.48573251d-02,
 462      1  1.70467655d-02, 1.93394396d-02, 2.17394852d-02, 2.42490773d-02,
 463      1  2.68684152d-02, 2.95957660d-02, 3.24275691d-02, 3.53586054d-02,
 464      1  3.83822285d-02, 4.14906520d-02, 4.46752791d-02, 4.79270545d-02,
 465      1  5.12368132d-02, 5.45956048d-02, 5.79949684d-02, 6.14271460d-02,
 466      1  6.48852271d-02, 6.83632267d-02, 7.18561029d-02, 7.53597274d-02,
 467      1  7.88708192d-02, 8.23868545d-02, 8.59059616d-02, 8.94268082d-02,
 468      1  9.29484864d-02, 9.64703968d-02, 9.99921344d-02, 1.03513375d-01,
 469      1  1.07033760d-01, 1.10552783d-01, 1.14069668d-01, 1.17583246d-01,
 470      1  1.21091827d-01, 1.24593066d-01, 1.28083828d-01, 1.31560049d-01,
 471      1  1.35016617d-01, 1.38447256d-01, 1.41844451d-01, 1.45199401d-01,
 472      1  1.48502033d-01, 1.51741065d-01, 1.54904135d-01, 1.57977973d-01,
 473      1  1.60948623d-01, 1.63801670d-01, 1.66522463d-01, 1.69096305d-01,
 474      1  1.71508595d-01, 1.73744902d-01, 1.75790974d-01, 1.77632682d-01,
 475      1  1.79255895d-01, 1.80646319d-01, 1.81789276d-01, 1.82669470d-01,
 476      1  1.83270725d-01, 1.83575716d-01, 1.83565712d-01, 1.83220322d-01,
 477      1  1.82517279d-01, 1.81432251d-01, 1.79938706d-01, 1.78007835d-01,
 478      1  1.75608540d-01, 1.72707519d-01, 1.69269456d-01, 1.65257378d-01,
 479      1  1.60633244d-01, 1.55358941d-01, 1.49398029d-01, 1.42718981d-01,
 480      1  1.35301474d-01, 1.27148627d-01, 1.18308730d-01, 1.08905085d-01,
 481      1  9.91559295d-02, 8.93515884d-02, 7.97824293d-02, 7.06663514d-02,
 482      1  6.21244732d-02, 5.41994827d-02, 4.68848207d-02, 4.01465202d-02,
 483      1  3.39357642d-02, 2.81954415d-02, 2.28635569d-02, 1.78750916d-02,
 484      1  1.31630892d-02, 8.65933391d-03, 4.29480447d-03/
 485       data u2 /
 486      1  7.17416019d-06, 1.70782645d-05, 3.31245126d-05, 6.01588363d-05,
 487      1  1.05339286d-04, 1.79174771d-04, 2.96719122d-04, 4.78862606d-04,
 488      1  7.53598916d-04, 1.15707860d-03, 1.73420412d-03, 2.53849668d-03,
 489      1  3.63099110d-03, 5.07800919d-03, 6.94782549d-03, 9.30645443d-03,
 490      1  1.22130079d-02, 1.57152366d-02, 1.98459102d-02, 2.46205841d-02,
 491      1  3.00370492d-02, 3.60764461d-02, 4.27057301d-02, 4.98809820d-02,
 492      1  5.75510102d-02, 6.56607602d-02, 7.41541974d-02, 8.29764928d-02,
 493      1  9.20754824d-02, 1.01402468d-01, 1.10912474d-01, 1.20564094d-01,
 494      1  1.30319039d-01, 1.40141489d-01, 1.49997326d-01, 1.59853293d-01,
 495      1  1.69676126d-01, 1.79431680d-01, 1.89084097d-01, 1.98595037d-01,
 496      1  2.07923034d-01, 2.17023055d-01, 2.25846345d-01, 2.34340694d-01,
 497      1  2.42451240d-01, 2.50121934d-01, 2.57297724d-01, 2.63927433d-01,
 498      1  2.69967170d-01, 2.75383917d-01, 2.80158840d-01, 2.84289739d-01,
 499      1  2.87792167d-01, 2.90698875d-01, 2.93057586d-01, 2.94927384d-01,
 500      1  2.96374262d-01, 2.97466488d-01, 2.98270390d-01, 2.98847025d-01,
 501      1  2.99249945d-01, 2.99524080d-01, 2.99705593d-01, 2.99822450d-01,
 502      1  2.99895431d-01, 2.99939301d-01, 2.99963931d-01, 2.99975129d-01,
 503      1  2.99974996d-01, 2.99961526d-01, 2.99927041d-01, 2.99854809d-01,
 504      1  2.99712769d-01, 2.99442742d-01, 2.98942676d-01, 2.98038511d-01,
 505      1  2.96441259d-01, 2.93684573d-01, 2.89040478d-01, 2.81421884d-01,
 506      1  2.69315148d-01, 2.50874185d-01, 2.24457680d-01, 1.89885662d-01,
 507      1  1.49894358d-01, 1.09927672d-01, 7.54041273d-02, 4.90259517d-02,
 508      1  3.06080023d-02, 1.85165524d-02, 1.09104125d-02, 6.27726960d-03,
 509      1  3.53002680d-03, 1.94049735d-03, 1.04218859d-03, 5.45964314d-04,
 510      1  2.77379128d-04, 1.33343739d-04, 5.32660444d-05/
 511       data u3 /
 512      1  1.86765383d-10, 1.96772458d-09, 1.19111389d-08, 5.54964761d-08,
 513      1  2.18340713d-07, 7.55899524d-07, 2.35604385d-06, 6.70801745d-06,
 514      1  1.76224112d-05, 4.30351929d-05, 9.82592148d-05, 2.10736217d-04,
 515      1  4.26209304d-04, 8.15657041d-04, 1.48160943d-03, 2.56186555d-03,
 516      1  4.22851247d-03, 6.68078970d-03, 1.01317466d-02, 1.47903961d-02,
 517      1  2.08424987d-02, 2.84336008d-02, 3.76573037d-02, 4.85502549d-02,
 518      1  6.10936693d-02, 7.52198901d-02, 9.08218891d-02, 1.07763660d-01,
 519      1  1.25889931d-01, 1.45034247d-01, 1.65025016d-01, 1.85689556d-01,
 520      1  2.06856371d-01, 2.28356037d-01, 2.50021072d-01, 2.71685149d-01,
 521      1  2.93181998d-01, 3.14344301d-01, 3.35002907d-01, 3.54986687d-01,
 522      1  3.74123404d-01, 3.92241969d-01, 4.09176451d-01, 4.24772089d-01,
 523      1  4.38893320d-01, 4.51433444d-01, 4.62324969d-01, 4.71549073d-01,
 524      1  4.79142163d-01, 4.85197409d-01, 4.89859810d-01, 4.93314543d-01,
 525      1  4.95770115d-01, 4.97439231d-01, 4.98520996d-01, 4.99187563d-01,
 526      1  4.99576941d-01, 4.99791928d-01, 4.99903753d-01, 4.99958343d-01,
 527      1  4.99983239d-01, 4.99993785d-01, 4.99997902d-01, 4.99999367d-01,
 528      1  4.99999835d-01, 4.99999965d-01, 4.99999995d-01, 5.00000000d-01,
 529      1  5.00000000d-01, 4.99999997d-01, 4.99999976d-01, 4.99999863d-01,
 530      1  4.99999315d-01, 4.99996914d-01, 4.99987300d-01, 4.99951740d-01,
 531      1  4.99829328d-01, 4.99435130d-01, 4.98245007d-01, 4.94883400d-01,
 532      1  4.86081966d-01, 4.65174923d-01, 4.21856650d-01, 3.47885738d-01,
 533      1  2.49649938d-01, 1.51648615d-01, 7.80173239d-02, 3.47983164d-02,
 534      1  1.38686441d-02, 5.05765688d-03, 1.71052539d-03, 5.38966324d-04,
 535      1  1.57923694d-04, 4.27352191d-05, 1.05512005d-05, 2.33068621d-06,
 536      1  4.45404604d-07, 6.88336884d-08, 7.23875975d-09/
 537       data zero/0.0d0/
 538 c
 539       abermx = zero
 540       do 10 k = 1,99
 541         ae1 = abs(y(1,k) - u1(k))
 542         ae2 = abs(y(2,k) - u2(k))
 543         ae3 = abs(y(3,k) - u3(k))
 544         abermx = max(abermx, ae1, ae2, ae3)
 545  10     continue
 546 c
 547       return
 548 c end of subroutine maxerr
 549       end
 550