trapezoid <- function(ftn, a, b, n = 100) {
  #trapezoidal numerical integration of ftn from a to b
  
  h <- (b - a)/n
  x.vec <- seq(a, b, by = h)
  f.vec <- sapply(x.vec, ftn)
  T <- h*(f.vec[1]/2 + sum(f.vec[2:n]) + f.vec[n+1]/2)
  return(T)
}