%
% *** function [ir,r] = newton(user_func,user_deriv,a,tol,itemax,iplot) ***
%
% Author: SJT
% Date: 1/19/2005
%
% Input
%   user_func  = an inline function defined by the user
%   user_deriv = an inline function (the derivative) defined by the user
%   a          = an estimate of the root
%   tol        = desired tolerance
%   itemax     = the maximum number of iterations allowed
%   iplot      = 1 for plotting
%
% Output
%   ir = 0 if Newton's method converges
%      = 1 if maximum number of iterations exceeded
%      = 2 if divide by zero
%   r  = an estimate of the root to the desired tolerance
%

function [ir,r] = newton(user_func,user_deriv,a,tol,itemax,iplot)

% Check input arguments
if nargin < 6, iplot = 0; end
if nargin < 5, itemax = 8; end
if nargin < 4, tol = 1E-6; end


% Set x1 and x2 to left-hand and right-hand ends of initial bracket
ir = 0;
da = 1;
cnt = 0;

while abs(da) > tol | abs(fa) > tol
    
    cnt = cnt + 1;
    fa = feval(user_func, a);
    df = feval(user_deriv,a);
    
    % Plot if desired
    if iplot == 1
        newton_plot(user_func,user_deriv,a,fa/df)
    end
    
    if abs(df) > 1E-13
        da = fa/df;
    else
        fprintf('ERROR: Divide by zero in function newton \n')
        ir = 2;
        return
    end
    a  = a - da;    
    fprintf('Iteration %3i, da = %12.6f, a = %12.6f, fa = %12.6f \n',...
             cnt, da, a, fa )

    % Store convergence data
    xstep(cnt) = da;
    
    % Check for maximum number of iterations
    if cnt >= itemax
        fprintf('ERROR: Maximum number of iterations exceeded \n')
        ir = 1;       
        break
    end
    
end

% Set best estimate of root
r = a;

fprintf('\nConvergence data for Newton method (quadratic) \n')
fprintf('xstep(i)/xstep(i-1)  xstep(i)/xstep(i-1)^2 \n') 
for i = 1:cnt-1
    xratio(i)  = xstep(i+1)/xstep(i);
    xratio2(i) = xstep(i+1)/xstep(i)^2;
    fprintf('%16.5e  %18.5e \n', xratio(i), xratio2(i))
end