program erfsim
c
c -- numerical integration by simpson's rule
c -- figure 11.4
c
      integer out
      real sum, upper, lower, tol
      external f
      common /inout/ out
      data lower/0.0/, upper/9.0/, tol/1.0e-5/
      data pi/ 3.141592/
c
      out = 6
      twopi = 2.0 / sqrt(pi)
10    write(out, 101)
      read(*, 102) upper
      if (upper .lt. 0.0) goto 99
      sum = 0.0
      if (upper .gt. 0.0)
     *   call simps(lower, upper, tol, sum, f)
      sum = sum * twopi
      write(out, 104) sum
      goto 10
99    stop
101   format(/' arg? ' )
102   format(e10.0)
104   format(' erf =', f10.5)
      end
      function f(x)
      f = exp(-x*x)
      return
      end
      subroutine simps(lower, upper, tol, sum, f)
c
c -- numerical integration by simpson's rule
c
      integer out, pieces, i, p2
      real x, delta, lower, upper, sum, tol
      real endsum, oddsum, sum1, evsum
      common /inout/ out
c
      pieces  = 2
      delta = (upper - lower) / pieces
      oddsum = f(lower + delta)
      evsum = 0.0
      endsum = f(lower) + f(upper)
      sum = (endsum + 4 * oddsum) * delta / 3.0
c     write(out, 101) pieces, sum
5     pieces = pieces * 2
       p2 = pieces / 2
       sum1 = sum
       delta = (upper - lower) / pieces
       evsum = evsum + oddsum
       oddsum = 0.0
       do 10 i = 1, p2
        x = lower + delta *(2 * i - 1)
        oddsum = oddsum + f(x)
10     continue
       sum = (endsum + 4.0 * oddsum + 2.0 * evsum)
     *    * delta / 3.0
c      write(out, 101) pieces, sum
      if (abs(sum - sum1) .gt. abs(tol * sum)) goto 5
      return
c 101   format(1x, i7, f9.5)
      end