function sol=newton(fun, x0, maxiter, tol)
% function sol = newton(fun, x0, maxiter, tol)
%
% find solution to system of n nonlinear equations
% with n unknowns of the form f(x) = 0 using Newton's method
% with numerical approximation of Jacobian (requires jacob.m)
%
% inputs:
%  fun - function whose root is being sought;
%       it is assumed that function accepts and returns column vector
%  x0 - initial pointas (column vector)
%  maxiter (optional) - maximum number of iterations (default 100)
%  tol (optional) - tolerance level (default eps); 
%		iteration stops if ||f(x_n)|| < tol
%
% outputs:
%  sol - numerical solution obtained
%
% (c) Grzegorz Klima 2004

if nargin < 4
  tol = eps;
end;

if nargin < 3
  maxiter = 100;
end;

x = x0;
n = 0;
f = feval(fun, x);

while ((f' * f >= tol) && (n < maxiter))
  j = jacob(fun, x);
  x = x - j \ f;
  f = feval(fun, x);
  n = n + 1;
end

if (f' * f >= tol)
  warning('newton did not converge...')
end

sol = x;