program simp2
c
c -- numerical integration by simpson's rule
c
      integer out
      real sum, upper, lower, tol
      external f, df
      common /inout/ out
      data lower/1.0/, upper/9.0/, tol/1.0e-5/
c
      out = 6
      write(out,*)' simpson''s rule integration'
      call simps(lower, upper, tol, sum, f, df)
      write(out, 104) sum
      stop

104   format(/' area =', f10.5/)
      end
      function f(x)
c
      f = 1.0 / x
      return
      end
      function df(x)
      df = -1.0 / (x * x)
      return
      end
c
      subroutine simps(lower, upper, tol, sum, f, df)
c
c -- numerical integration by simpson's rule
c -- with end correction
c
      integer out, pieces, i, p2
      real x, delta, lower, upper, sum, tol
      real endsum, oddsum, sum1, evsum, endcor
      common /inout/ out
c
      pieces  = 2
      delta = (upper - lower) / pieces
      oddsum = f(lower + delta)
      evsum = 0.0
      endsum = f(lower) + f(upper)
      endcor = df(lower) - df(upper)
      sum = (endsum + 4 * oddsum) * delta / 3.0
      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 = (7 * endsum + 16 * oddsum + 14 * evsum
     *    + endcor * delta) * delta / 15.0
       write(out, 101) pieces, sum
      if (abs(sum - sum1) .gt. abs(tol * sum)) goto 5
      return
101   format(1x, i7, f9.5)
      end