function x = syszero(fsys,fjac,x0,xtol,rtol,maxit)
% x = syszero(fsys,fjac,x0,xtol,rtol,maxit)
%
% Naive Newton method for nonlinear systems.
%
% Inputs:
%   fsys   String, name of function that returns residual vector
%   fjac   string, name of function that returns the Jacobian
%   x0     starting guess vector
%   xtol   stopping tolerance for change in x
%   rtol   stopping tolerance for residual
%   maxit  maximum # of iterations permitted
%   
% Classroom example:
% >> x = syszero('F3d','dF3d',[0;2;2],1e-14,1e-14,10)
% 
x = x0(:);                         % Column vector
disp('')                           % Header for output table
disp('  m    dx     || F(x) ||')
disp('-------------------------')
r = feval(fsys,x);  % Initial residual
s = 1 + xtol;       % Guarantee at least one iteration :-)
m = 0;              % # iterations completed
while ( norm(r)>rtol & norm(s)>xtol & m<maxit )
  J = feval(fjac,x);
  s = -J\r;
  x = x+s;
  r = feval(fsys,x);
  m = m + 1;
  disp(sprintf('%3d  %8.2e  %8.2e',m,norm(s),norm(r)))
end