%
function [x] = newton(f, fp, x0, tolerance, maxiter, printflag)
% newton  newton's method
% f = function
% fp it's derivative
% x0 = initial guess
% tolerance = maximum error desired
% maxiter = maximum iterations
% printflag = true, display intermediate results
% printflag = false, display only the final result

   x = x0;
   
   if printflag
      fprintf('%3s%24s%24s\n', 'n','x','err');
   end
   
   n = 1;
   while n <=  maxiter
      fpx = fp(x);
      
      if fpx == 0
         error('division by zero');
      end
      
      dx = f(x)/fpx;
      x = x - dx;
      err = abs(dx);
      
      if printflag
         fprintf('%3d%24.15e%24.15e\n', n,x,err);
      end
      
      if err < tolerance
         break;
      end

   n = n + 1;
   end
   
   if n > maxiter
      fprintf('Failure to converge in %d iterations\n', maxiter);
   end
   % x is automatically returned here
end