subroutine integrate(f,a,b,n,val)
  real, intent(in) :: a, b
  real, intent(out):: val
  real, external   :: f     ! function
  integer          :: i, n
  real             :: h, sum
  h = (b-a)/n
  ! position 1
  sum = 0.5*f(a)
  ! position 2
  sum = sum + 0.5*f(b)
  do i = 1, n-1
    ! position 3
    sum = sum + f(a+i*h)
  end do
  ! position 4
  val = sum*h;
end subroutine integrate