function x=Newton(f,x0,tol,maxiter,fp);
% NEWTON solves a nonlinear system of equations   
%   x=Newton(f,x0,tol,maxiter,fp); solves the nonlinear system of
%   equations f(x)=0 using Newtons methods, starting with the initial
%   guess x0 up to a tolerance tol, doing at most maxit iterations. An
%   analytical Jacobian can be given in the parameter fp
    
numJ=nargin<5;
if nargin<4, maxit=100; end;
if nargin<3, tol=1e6; end;

x=x0; i=0; dx=ones(size(x0));   
while norm(f(x))>tol & norm(dx)>tol & i<maxiter
  if numJ, J=NumericalJacobian(f,x); else J=fp(x); end;
  dx=-J\f(x); x=x+dx; i=i+1;
end
if i>=maxiter 
  error('Newton did not converge: maximum number of iterations exceeded'); 
end;
