function x = NewtonD(F,JF,x,tol,maxIt)
% D-dimensional Newton'n method to
% solve the system F(x) = 0
% F(x)=[F1(x); F2(x); ... ; Fn(x)]
% x = [ x1; x2; ... ; xn]
% JF is the Jacobian matrix of F at x
% x is the initial guess
% tol is the tolerance
% maxIt is the maximum number of iterations
% Usage: x = NewtonD('fcnEg5p15','fcnJEg5p15',[1;2],1e-6,20)
% K. Ming Leung, 01/03/03

disp([0 x']);
iter=1; dif=10;
while (dif > tol & iter <=maxIt)
    dx = -feval(JF,x) \ feval(F,x);
    dif = norm(dx);
    x = x + dx;
    disp([ iter x' dif]);
    iter = iter + 1;
end;
if iter > maxIt
    disp('Maximum iteration exceeded without convergence!');
else
    disp('Newton method has converged.');
end

function f = fcnEg5p15(x)
f = [(x(1)+2*x(2)-2); (x(1)^2+4*x(2)^2-4)];

function Jf = fcnJEg5p15(x)
Jf = [1 2; 2*x(1) 8*x(2)];

function f = orbitCross(x)
f = [ (3*x(1)^2 + 4*x(2)^2 -3)
(x(1)^2 + x(2)^2 - sqrt(3)/2)];

function df = orbitCrossJ(x)
df = [6*x(1) 8*x(2); 2*x(1) 2*x(2)];