program csr
  ! Purpose: to calculate the integral of f(x) over the interval [a,b],
  ! using the composite Simpson's rule with s panels (n=2*s subintervals
  ! or 2*s+1 points). The 2*s+1 points are numbered from 0 to 2*s.
  !    Input:  a,b - endpoints of integration interval
  !    Output: s, simp - number of panels, approximation to integral
  !    Subprogram required: f(x) - returns value of integrand at x.

implicit none
integer :: s, n, i, k
real :: a, b, f, h, simp, sum1, sum2, sum3, t
character(len=1) :: check
 
problem_loop: do

   ! Read in and echo a,b:
   write(*,*) 'Input endpoints a,b:'
   read*, a,b
   print*, 'Endpoints a,b:', a,b
   ! Calculate contribution from end-points:
   sum1 = f(a) + f(b)   
   ! Print header:
   write(*,*) 'Number of panels  Compound Simpson''s rule'
   s = 1
   panel_loop: do k = 1,10
      ! Calculate number of subintervals n:
      n = 2*s
      ! Compute subinterval size:
      h = (b - a)/n 
      ! Calculate contribution from even points:
      sum2 = 0.
      if(n/=2) then
         do i = 2,n-2,2
            t = a + i*h
            sum2 = sum2 + f(t)
         enddo
      endif 
      ! Calculate contribution from odd points:
      sum3 = 0.
      do i = 1,n-1,2
         t = a + i*h
         sum3 = sum3 + f(t)
      enddo
      ! Calculate the integral:
      simp = h*( sum1 + 2.*sum2 + 4.*sum3)/3.

      write (*,10) s, simp
10    format(6x,i5,10x,e13.6)
      ! Double number of panels:
      s = 2*s
   enddo panel_loop

   write(*,20)
20 format(//,'Do you wish to continue ? (y/n) :')
   read (*,30) check
30 format(a)
   if(check/='Y'.and.check/='y') stop
enddo problem_loop

end program csr	
!=======================================================
function f(x) 
! Integrand.
implicit none
real :: f, x
f = 4.0/(1.0+x*x)
return
end function f