Program main
!=============================================
! Integration of a function using Simpson rule 
!=============================================
implicit none
double precision f, a, b, integral
integer n, i
double precision, parameter:: pi = 3.1415926
external f

a = 0.0
b = pi

n = 2

write(*,100)

do i=1,16
   call simpson(f,a,b,integral,n)
   write (*,101) n, integral
   n = n*2
end do

100   format('     nint   Simpson')
101   format(i9,1pe15.6)
end


  Function f(x)
!----------------------------------------
! Function for integration
!----------------------------------------
implicit none
double precision f, x
 f = sin(x)
! f = x*cos(10.0*x**2)/(x**2 + 1.0)
return
end

 Subroutine simpson(f,a,b,integral,n)
!==========================================================
! Integration of f(x) on [a,b]
! Method: Simpson rule for n intervals  
! written by: Alex Godunov (October 2009)
!----------------------------------------------------------
! IN:
! f   - Function to integrate (supplied by a user)
! a	  - Lower limit of integration
! b	  - Upper limit of integration
! n   - number of intervals
! OUT:
! integral - Result of integration
!==========================================================
implicit none
double precision f, a, b, integral,s
double precision h, x
integer nint
integer n, i

! if n is odd we add +1 to make it even
if((n/2)*2.ne.n) n=n+1

! loop over n (number of intervals)
s = 0.0
h = (b-a)/dfloat(n)
do i=2, n-2, 2
   x   = a+dfloat(i)*h
   s = s + 2.0*f(x) + 4.0*f(x+h)
end do
integral = (s + f(a) + f(b) + 4.0*f(a+h))*h/3.0
return
end subroutine simpson