function [x,n] = newton(f,fd,x0,atol)
%
% This function implements Newton's method
%      x_k+1=x_k - f(x_k)/fd(x_k)
% written by Jessica Bosch
%
% x0 is the initial guess, f the function for which we want to find the
% root, and fd the derivative of f.
% Upon success, this function returns in x a value x_k such that
%      | x_k - x_k-1 | < atol,
% and in n the number of iterations required.

    
  % initialization    
    x_old=x0+1;                        % x_old will be x_k-1
    x_new=x0;                          % x_new will be x_k
    fx = feval(f,x_new);               % f(x_0)
    fxd = feval(fd,x_new);             % f'(x_0)
    n=0;                               % iteration counter


  % Newton's method
    while (abs(x_new-x_old)>=atol)

        x_old = x_new;              % update previous iterate      
        x_new = x_new - fx / fxd;   % compute next iterate
        
        fx = feval(f,x_new);        % f(x_k+1)
        fxd = feval(fd,x_new);      % f'(x_k+1)
        n=n+1;                      % update iteration counter

        % break if it takes too long
        if (n>=1e5)
           break;
        end
    end

    x=x_new;                        % approximate root
    
end
