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