find.mle <- function(f0,f1,f2,initial.estimate,tolerance=0.01,show.plots=TRUE,range=c(0,1),pause=TRUE) {
  if ( show.plots && pause ) {
    old.ask <- par()$ask
    par(ask=TRUE)
  }
  estimate <- initial.estimate
  iter <- 0
  while ( abs(f1(estimate)) > tolerance ) {
    b <- f1(estimate)
    c <- f2(estimate)
    new.estimate <- estimate - b/c
    if ( show.plots ) {
      a <- f0(estimate)
      x <- seq(range[1],range[2],by=0.01)
      ll <- f0(x)
      plot(x,ll,ylab="f(x)",type="l")
      abline(v=estimate)
      abline(v=new.estimate,col="red")
      taylor.approx <- a + b*(x-estimate) + 0.5*c*(x-estimate)^2
      lines(x,taylor.approx,col="green")
    }
    estimate <- new.estimate
    iter <- iter + 1
  }
  cat("It took",iter,"iterations to converge.\n")
  if ( show.plots && pause ) {
    par(ask=old.ask)
  }
  estimate
}