function r = newton( fun, funder, xinit, tol, nmax, verbose)

% newton  Uses Newton's method to find a root 
%         of the equation  fun(x) = 0 
%         starting at xinit as an initial approximation
%
% Input:  fun = the function whose zero we are searching
%         funder = the derivative of the function fun
%         xinit = the value from which we start the search 
%         tol = x tolerance
%         nmax = maximum number of iterations
%         verbose = 1 in case all iterations are needed, 
%                   0 otherwise
%
% Output: r = the root of fun(x) = 0 with the desired tolerance

if tol < 3*eps       %  The tolerance is too small
  warning(sprintf('the tolerance is too small\n'));
  return;
end

if verbose == 1         %  Printing the headings
  fprintf(['\niters         x                 x-xold' ...
	   '    log_10(|x-xold|)\n\n']);
end


%  The main part of the program starts here

xold = xinit;

for i = 1:nmax
  fx = feval(fun,xold);       %  fx = fun(xold)
  dfdx = feval(funder,xold);  %  dfdx = fun'(xold)

  x = xold - fx/dfdx;         %  Main step in Newton's method

  if verbose == 1             %  Printing all iterations
    fprintf('%3d  %18.14f  %18.14f  %10.5f\n', ...
		        i, x, x-xold,log(abs(x-xold))/log(10.) );
  end
      
  if abs(x-xold) < tol
    r = x;    return;      %  We found the root
  end

  xold = x;
	  
end

warning(sprintf('root not found within tol after %d steps\n',i));