/* -*-ePiX-*- */
#include "epix.h"
using namespace ePiX;

double f(double t)
{
  return 1 + t*(-3 + t*t);
}

double df(double t)
{
  return Deriv(f)(t).x2();
}

int main()
{
  picture(P(1.5,0), P(2,3), "2.5 x 2.5in");

  begin();
  h_axis(2);

  double x0(2),  x1(x0-f(x0)/df(x0));
  label(P(x0,0), P(0,-4), "$x_i$", b);

  for (int i=0; i < 3; ++i)
    {
      dashed();
      line(P(x0,0), P(x0,f(x0)));
      solid();
      line(P(x0,f(x0)), P(x0-f(x0)/df(x0),0));

      x0 -= f(x0)/df(x0); // Newton's method
    }

  bold();
  plot(f, xmin(), xmax()+0.05, 60);

  label(P(1.75, f(1.75)), P(-2,2), "$y=f(x)$", tl);
  label(P(x1,0), P(0,-4), "$x_{i+1}$", b);
  label(P(1.75, df(2)*(1.75-x1)), P(0,-4), "Slope $= f'(x_i)$", br);

  end();
}