function ShowNewton()
% Illustrates the finite-difference Newton idea for systems of nonlinear
% equations.
x0 = [1;1;1];
tol = .000000000001;
itMax = 20;
xstar = FDNewton(@MyF,@MyJ,x0,tol,itMax);

function xStar = FDNewton(F,J,x0,tol,itMax)
% Finds a zero of the function F:Rn ---> Rn with starting vector x0.
% Uses the Finite Difference Newton Method.
% Terminates as soon as norm(F(xStar))<=tol or after itMax steps, whichever
% comes first.
xc = x0;
Fc = F(xc);
its = 0;
while norm(Fc)>tol && its<itMax
Jc = J(F,xc,Fc);
sc = -(Jc\Fc);
xc = xc + sc;
Fc = F(xc);
its = its+1;
fprintf('%3d %10.3e %20.15f %20.15f %20.15f\n',its,norm(Fc),xc)
end
xStar = xc;

function Fc = MyF(xc)
% xc is a column 3-vector
% This function is zeroed by a point that is on three given elliposids
A1 = [6 -7 2; -7 10 -3; 2 -3 1];
A2 = [5 -1 1; -1 10 -1; 1 -1 1];
A3 = [9 4 -2; 4 2 -1; -2 -1 1];
v = [2;1;-1];
Fc = [(xc-v)'*A1*(xc-v)-1;...
(xc-v)'*A2*(xc-v)-1;...
(xc-v)'*A3*(xc-v)-1];

function Jc = MyJ(F,xc,Fc)
% The Jacobian of the function F at xc. Fc = F(xc)
n = length(xc);
Jc = zeros(n,n);
for k=1:n
% The divided difference parameter...
delta = sqrt(eps)*xc(k);
ek = [zeros(k-1,1);1;zeros(n-k,1)];
% Approximate the partial of F(x) with respect to x(k)
Jc(:,k) = (F(xc+delta*ek)-Fc)/delta;
end