function trapezoid(a,b,n) result(s)
    real,    intent(in) :: a
    real,    intent(in) :: b
    integer, intent(in), optional :: n
    real :: s
    real :: h
    integer :: i, div = 10
    if (present(n)) div = n
    h = (b-a) / div
    s = (f(a)+f(b)) / 2.0
    do i = 1, n-1
      s = s + f(a + i*h)
    end do
    s = h*s
  end function trapezoid