simpson <- function(ftn, a, b, tol = 1e-8, print_n = FALSE) {
  #Simpson numerical integration of ftn from a to b with tolerance
  
  #initialize
  n <- 4
  h <- (b-a)/4
  fx <- sapply(seq(a, b, by =h), ftn)
  S <- sum(fx*c(1, 4, 2, 4, 1))*h/3
  S.diff <- tol + 1
  
  #increase n
  while(S.diff > tol) {
    S.old <- S
    n <- 2*n
    h <- h/2
    fx[seq(1, n+1, by = 2)] <- fx
    fx[seq(2, n, by = 2)] <- sapply(seq(a+h, b-h, by = 2*h), ftn)
    S <- h/3*(fx[1] + 4*sum(fx[seq(2, n, by = 2)]) + 
                2*sum(fx[seq(3, n-1, by = 2)]) + fx[n+1])
    S.diff <- abs(S - S.old)
  }
  if (print_n) cat("andaze baze ha", n, "\n")
  return(S)
}
  
  