function x = newton(func,x0,param,crit,maxit);
%NEWTON   Finds the fixed point of a system of equations.
%         x=newton(func,x0,param,crit,maxit)  uses the Newton-Raphson
%         algorithm to find the zeros of the following system:
%
%                        f1(z1,z2,...zn)=0
%                        f2(z1,z2,...zn)=0                (*)
%                          :           :
%                        fn(z1,z2,...zn)=0
%
%         where x=[z1,z2,...,zn]' are those values of zi that solve (*).
%         To use newton, you must specify "func" which is a string giving
%         the name of an m-file.  Do not use "f.m","fj.m","x.m","x0.m".
%         This m-file takes inputs x and param and returns f and fj where:
%
%           f = [f1(x,param),f2(x,param),...fn(x,param)]'  
%           J = [df1(x,param)/dz1 df1(x,param)/dz2 ... df1(x,param)/dzn;
%                df2(x,param)/dz1 df2(x,param)/dz2 ... df2(x,param)/dzn;
%                       :                 :                 :
%                dfn(x,param)/dz1 dfn(x,param)/dz2 ... dfn(x,param)/dzn]
%
%         The arguments x0, crit, and maxit are the initial guess for the
%         fixed point, the convergence criterion, and the iteration limit.

%         Ellen R. McGrattan, 11-13-87
%
if nargin < 4; crit=1e-8; maxit=100; end;
n=length(x0);
for i=1:maxit;
  eval(['[f,fj]=',func,'(x0,param);']);
  x=x0-fj\f;
  if max( abs(x-x0)./(1+abs(x0)) )<crit; 
    break; 
  end;
  x0=x;
end;
if i>=maxit; 
  sprintf('WARNING: iteration limit of %g exceeded for NEWTON',maxit)
end;
