function [x, hist] = newtsol(x, f, tola, tolr,jdiff) 
%NEWTSOL   Newton-Armijo code for scalar equations.
%
% This code terminates on small relative-absolute errors
%
%         [X, HIST] = NEWTSOL(X, F, TOLA, TOLR, JDIFF) 
%
% Inputs: X = initial iterate
%         F = function
%         TOLA = absolute error tolerance
%         TOLR = relative error tolerance
%         JDIFF = 0, analytic derivative 
%                    syntax: [f,fp]=f(x)
%         JDIFF = 1, forward difference derivative 
%
% Output: x=approximate result
%      HIST = array of iteration history, useful for tables and plots
%            The four columns are iteration number, residual, 
%            number of step size reductions done in the line
%            search, and the Newton iterate itself
%
%
%             If you leave out the hist array by calling newtsol as
%             Z = NEWTSOL(X, F, TOLA, TOLR) 
%                then storage for the history array is not allocated and
%                the iteration history is not stored.
%
if nargin == 4
    jdiff = 1;
end
maxit=100;    % maximum iterations
maxitls=20;   % maximum iterations inside the line search
alpha=1.d-4; 
small=1.d-8;
%
% Initialize
%
itc=0;
fc=feval(f,x);
tol=tolr*abs(fc)+tola;
%
% Store iteration history?
%
if nargout == 2
    histmp=zeros(maxit,4);
    histmp(itc+1,:)=[itc, abs(fc), 0, x];
end
%
%    Main Newton loop
%
while(abs(fc) > tol)
%
%    Iteration becoming unbounded?
%
    if(abs(x) > 1.d9)
        disp(' iterate too large');
        if nargout == 2
            histmp(itc+1,:)=[itc, abs(fc), 0, x];
            hist=histmp(1:itc+1,:);
        end
        return
    end
%
    lambda=1;
    iarm=0;           % line search iteration counter
    if jdiff==1
        fp=fordiff(f,x,fc);% use a numerical derivative
    else
        [fv,fp]=feval(f,x);
    end
%
%    Check for small relative derivatives.
%
    if(abs(fp) <= small*abs(fc))
        disp(' small derivative error');
        return;
    end
    s=-fc/fp;
    xt=x+lambda*s;
    ft=feval(f,xt);
%
%    Main line search loop. 
%
    while(abs(ft) >= (1 - alpha*lambda)*abs(fc) + 1.d-12)
        lambda=lambda/2;
        xt=x+lambda*s;
        ft=feval(f,xt);
        iarm=iarm+1;
%
%    Are you spending too much time in the line search?    
%
        if(iarm > maxitls)
            disp(' line search failure');
            if nargout == 2
                histmp(itc+1,:)=[itc, abs(fc), 0, x];
                hist=histmp(1:itc+1,:);
            end
            return
        end
    end
%
%    Step accepted, continue the Newton iteration.
%
    x=xt;
    itc=itc+1;
%
%    Too many iterations?
%
    if(itc > maxit)
        disp(' maxit reached');
        if nargout == 2
            histmp(itc+1,:)=[itc, abs(fc), 0, x];
            hist=histmp(1:itc+1,:);
        end
        return
    end
%
    fc=feval(f,x);
    if nargout == 2
       histmp(itc+1,:)=[itc,abs(fc),iarm,x];
    end
end
%
%   Fix up the history array.
%
if nargout == 2
    hist=histmp(1:itc+1,:);
end
%
% A simple forward difference in MATLAB.
%
function fp = fordiff(f,x,fc)
h=1.d-7;
%
%       Notice how f is evaluated using the feval command
%
fright=feval(f,x+h);
fp=(fright-fc)/h;