real function trapezoidal(a,b,n)
  real :: a, b, f, sum, h
  external :: f
  integer :: i, n
  h = (b-a)/n
  sum = 0.5*(f(a) + f(b))
  do i = 1, n-1
     sum = sum + f(a+i*h)
  end do
  trapezoidal = h*sum
end function trapezoidal