function[x] = newton(F, x0, iterations)
% Function NEWTON(F, x0, iterations)
% Estimate zeros of F by Newton's method.
%  F should be a column vector n high of
% real-valued functions of n variables.
% These variables must be named x1 ... xn.
%  x0 is an initial guess, a real column n-vector.
%  iterations is the number of iterations to compute.
% EXAMPLE:
%  syms x1 x2; newton([cos(x1)-x1; sin(x2)], [.1; 3.0], 3)

%
% Calculate the Jacobian matrix (derivative) symbolically.

vars = '[';

% construct an array of symbolic variable names for
% Jacobian 

for i = 1:length(F)
    iS = num2str(i);
    vars = [vars 'x' iS ' '];
    eval(['x' iS ' = sym(''x' iS ''');']);
end
vars = [vars ']' ];
eval(['vars= ' vars ';']);
J = jacobian(F, vars);

% call Maple routine
% Iterate.

x = x0;
for i = 1:iterations
    JJ = double(subs(J, vars, x.'));
    FF = double(subs(F, vars, x.'));
    x = x - inv(JJ) * FF
end
