newton <- function(f3, x0, tol = 1e-9, n.max = 100) {
  #newton algorithm in multivariate optimization
  
  x <- x0
  f3.x <- f3(x)
  n <- 0
  while ((max(abs(f3.x[[2]])) > tol) & (n < n.max)) {
    x <- x - solve(f3.x[[3]], f3.x[[2]])
    f3.x <- f3(x)
    n <- n + 1
  }
  if (n == n.max) {
    cat("newton failed to converge\n")
  } else {
    return(x)
  }
}
