% WARNING!  This is for f(x) = (x+1)(x-1)^2.  If you're going to
% change the function, and its derivative in _two_ places!
%
% [root,value,steps,A] = newton(x_0,tol) reads in a first guess
% x_0 and a desired tolerance tol.  it
% takes steps until it gets |f(root)| < tol.  It then returns the root,
% the value of the function at the root, and the number of steps it
% took to find it

% A is a diagnostic matrix.  At the moment, it's 100x3 and contains
% A(i,1) = i, A(i,2) = x_i, A(i,3) = f(x_i)

function [root,value,steps,A] = newton(x_0,tol)

f = 2*tol;
steps = 0;

A = zeros(100,3);

while ( abs(f) > tol) & (steps < 100) 

% first evaluate f and f_x at the present point (really it's x_i, but we
% call it x_0
% this is where you'll have to change f and f_x!!!
    f = (x_0+1)*(x_0-1)^2;
    f_x = (x_0-1)^2 + 2*(x_0+1)*(x_0-1);

    steps = steps + 1;
    A(steps,1) = steps;
    A(steps,2) = x_0;
    A(steps,3) = f;	  

% now use this to define the next point (really it's x_(i+1), but we
% call it x_0   x_(n+1) = x_n - f(x_n)/f_x(x_n):
    x_0 = x_0 - f/f_x;

% print to the screen the step number and the function size
    steps, f

end

root = x_0;
value = f;	
