%1.2 Subroutine newtons
function [x, fx, xx] = newtons(f,x0,TolX, MaxIter,varargin)
% newtons.m: Newton Raphson root zeros for nonlinear system of equations
% input:
% f = 1^st order function vector
% x0 = the initial guess of the solution
% TolX = convergence criteria
% MaxIter = maximum number of iterations
% output:
% x = the point at which the algorithm has reasched
% fx = the value of the function at x
% xx = the history of the function
h = 1e-5; TolFun = 1e-5; EPS = 1e-6; fx = feval(f,x0,varargin{:});
Nf = length(fx); Nx = length(x0);
if Nf~=Nx, error('Incompatible dimensions of f and x0'); end
if nargin <4, MaxIter = 100; end 
if nargin < 3, TOLx = EPS; end
xx(1,:) = x0(:)'; % Initialize the solution as a row vector

for k=1:MaxIter
dx = -jacobi(f,xx(k,:),h,varargin{:})\fx(:); % Newton method [dfdx]^-1*fx
xx(k+1,:) = xx(k,:) + dx';
fx = feval(f,xx(k+1,:),varargin{:});
if norm(fx) < TolFun | norm(dx) < TolX, break; end
end
x = xx(k+1,:);
if k == MaxIter, fprintf('The best in %d iterations \n', MaxIter),
end