program ex032
c This program integrates function f1(=pi*sin(pi*x)) from 0 to 1
c by using the Simpson formula
      implicit real*8 (a-h,o-z)
      external f1
      a=0.0d0
      b=1.0d0
      exa=2.0d0
          do n=10,100,10
c     call daikeif(f1,a,b,n,sum)
      call simpsonf(f1,a,b,n,sum)
      dx=(b-a)/n
      err=abs(sum-exa)
      write( *,800) n, dx, sum, exa, err, err/dx/dx/dx/dx
      write(16,800) n, dx, sum, exa, err, err/dx/dx/dx/dx
          enddo
 800  format(1x,i4,5(1x,e11.4))
      end
c **********************************************************************
      subroutine simpsonf(FUNC,a,b,n,sum)
c **********************************************************************
c
c     This subroutine calculate the integral
c       sum = integral_a^b FUNC(x) dx
c     by using the Simpson formula.
c
      implicit real*8 (a-h,o-z)
      external FUNC
c
      sum = 0.0D0
      dx = (b - a)/n
      do i = 1,n
         x1 = a + (i - 1.0D0)*dx
         x2 = a + (i - 0.5D0)*dx
         x3 = a +  i         *dx
         sum = sum + (FUNC(x1)+4*FUNC(x2)+FUNC(x3))/6
      end do
      sum = sum*dx
c
      end
c **********************************************************************
      function f1(x)
c **********************************************************************
c
c     function f1(x)=pi*sin(pi*x)
c
c **********************************************************************
      implicit real*8 (a-h,o-z)
      pi=4.0d0*atan(1.0d0)
c
      f1=pi*sin(pi*x)
c
      end