function [x,val,iter,flag]=newton(fun,dfun,xstart,tol,maxiter)
%[x,val,iter,status]=newton(fun,xstart,tol,maxiter)
% NEWTON: zero finding with Newton's method
%
% Input:      fun: the function (function handle or name)
%            dfun: a function which computes the derivative of fun
%          xstart: initial guess for the zero
%             tol: tolerance for abortion:     |fun| < tol -> success
%                                             |dfun| < tol -> critical point         
%                                         |dfun\fun| < tol -> stagnation
%                  (default 1e-8)                
%         maxiter: maximum iteration cout (default 100)
%
% Output:       x: the last approximation for the zero 
%             val: the value of fun at x
%            iter: the number of iterations done 
%            flag: 1: |fun(x)| < tol
%                 -1: maxiter reached
%                 -2: |dfun(x)| < tol 
%                 -3: |dfun(x)\fun(x)| < tol 

% Author:  Joerg Hoerner

% input check and defaults
if nargin < 3 
  error('newton needs at least three input arguments');
end

if ~exist('tol','var') | isempty(tol)
  tol=1e-8;
end
if ~exist('maxiter','var') | isempty(maxiter)
  maxiter=100;
end

% initialisation
x=xstart;
val=feval(fun,x);
iter=0;
flag=0;
	
% iteration
while (norm(val) > tol & iter < maxiter)
  iter=iter+1;
  df=feval(dfun,x);   % derivative
  if norm(df) < tol,flag=-2;return,end    % critical point
  delta=-df\val;       
  if norm(delta) < tol,flag=-3;return,end % stagnation
  x=x+delta;          % update
  val=feval(fun,x);     
end

% set status
if iter==maxiter
  flag=-1;
else
  flag=1;
end



