function [x] = newton(f, df, x0, tol) 
%
% Computed the root of f(x)=0 by using Newton's method
% Input:
%    f: nonlinear function
%   df: the derivative of the function f
%   x0: initial value
%  tol: stop tolerance
% Output:
%    x: approximate solution
%
max_it = 50;
it     = 1; 
err    = 1; 
fprintf(' k              x_k                |x_k-x_{k-1}| \n');
fprintf('-------------------------------------------------\n');
fprintf('%2.0f    %24.16e   \n',it-1, x0);
while ( it <= max_it && err >= tol )
    x   = x0 - feval(f,x0) / feval(df,x0);
    err = abs(x - x0);
    it  = it + 1; 
    x0  = x;
    fprintf('%2.0f    %24.16e   %14.4e \n',it-1, x, err);
end