function  [x,xhist,fhist] = newton(fname,x0,S1,P1,S2,P2,S3,P3,S4,P4)

% NEWTON  Find a zero of a function whose derivative is available.
%    X = NEWTON('FNAME',X0) finds a zero of a function whose name
%    is given in the string 'FNAME', using Newton's method with
%    starting guess X0. The function m-file FNAME takes one input X
%    and returns two outputs, the function value F and the function 
%    derivative (jacobian)  J at X.
%    X and F may be complex and may be vector-valued.
%    NEWTON('FNAME',X0,Param1,Val1,..ParamN,ValN) allows specification of
%    additional parameters. The parameters can be specified
%    in any order or may be omitted. The parameters are:
%      'XTol' : Iterations stop when the 2-norm of the difference between
%         successive X values is less than this value. Default = 1e-10.
%      'FTol' : Iterations stop when the 2-norm of the function value
%         is less than this value. Default = 1e-10.
%      'MaxSteps' : Maximum allowed number of iterations. Default = 10.
%      'FPar' : Constant parameters used by function; function call
%          is then FNAME(X,Fpar).
%    [X,XHIST,FHIST] = NEWTON(...) returns the X and F values for 
%    successive iterations, one column per iteration.
%    For example, NEWTON('myfunc',3,'FTol',1e-6) finds a zero near X0=3.

% Ref: Kincaid & Cheney, Numerical Analysis, Brooks-Cole, 1991, p.64.
% Author: Robert Piche, Tampere University of Technology, Finland
% Last updated: Nov. 26, 1993

if nargin==0, help newton, return, end
%===============================================================
% Default values
%===============================================================
XTol = 1e-10;
FTol = 1e-10;
MaxSteps = 10;
FParExists = 0;

%===============================================================
% Set user-defined parameters
%===============================================================
parms = ['XTol    '
         'FTol    '
         'MaxSteps'
         'FPar    '];
[n1,n2]=size(parms);
nn = nargin;
ni = 1;
while nn > 2
  Sin = eval(['S' int2str(ni)]);
  S = [Sin blanks(n2-length(Sin))];
  if any(all(S(ones(n1,1),:)'==parms'))
    eval([ S '= P' int2str(ni) ';']);
    if S==parms(4,:), FParExists = 1; end
  else
    disp(['WARNING: unknown parameter ''' Sin ''' for NEWTON.'])
  end
  ni = ni+1;
  nn = nn - 2;
end

%===============================================================
% Newton's method
%===============================================================
if FParExists
  [v,J] = feval(fname,x0,FPar); 
else 
  [v,J] = feval(fname,x0); 
end
xhist = x0; 
fhist = v; 
x = x0;
if norm(v) < FTol, return, end
for n=1:MaxSteps
  dx = J\v;
  x  = x - dx;
  if FParExists
    [v,J] = feval(fname,x,FPar); 
  else 
    [v,J] = feval(fname,x); 
  end
  xhist = [xhist x ]; 
  fhist = [fhist v ];
  if ( norm(dx) < XTol )|( norm(v) < FTol ), return, end
end

if (XTol>0)|(FTol>0)
  disp('WARNING: NEWTON reached max. number of iterations without converging')
end