function [x, hist] = secant(x, f, tola, tolr) 
%
%SECANT  This a secant code with no line search. 
%
% This code terminates on small relative-absolute errors
%
% function [X, HIST] = SECANT(X, F, TOLA, TOLR) 
%
% Inputs: X = initial iterate
%         F = function
%         TOLA = absolute error tolerance
%         TOLR = relative error tolerance
%
% Output: X = approximate result
%      HIST = array of iteration history, useful for tables and plots
%                The two columns are iteration number and residual.
%
%             If you leave the hist array  out by calling secant as
%             Z = SECANT(X, F, TOLA, TOLR) 
%             the storage for hist is not allocated and
%             the iteration history is not stored.
%
% Limit iterations to 100.
%
maxit=100;    
if nargout == 2
    histmp=zeros(maxit,2);
end
%
% Initialize
%
itc=0;
fc=feval(f,x);
%
% A secant code needs two iterations and the
% function values at each. I'll start it off with the
% old value xm = .99 x, if x ~= 0, or xm=.001, if x = 0,
% where x is the initial iterate.
%
if x~=0
    xm=x*.99;
else
    xm=.001;
end
fm=feval(f,xm);
%
tol=tolr*abs(fc)+tola;
%
%    Store iteration history?
%
if nargout == 2
    histmp(itc+1,:)=[itc, fc];
end
%
%    Main Secant loop.
%
while(abs(fc) > tol)
%
%    Iteration becoming unbounded?
%
    if(abs(x) > 1.d9)
        disp(' iterate too large');
        if nargout == 2
            histmp(itc+1,:)=[itc, fc];
            hist=histmp(1:itc+1,:);
        end
        return
    end
%
% Use a secant approximate derivative.
%
    del=x-xm;
    if del ~= 0
        fp=(fc-fm)/del;
        s=-fc/fp;
    else
        s=0;
    end
%
% Continue the secant iteration. Make sure you update xm and fm.
%
    fm=fc; xm=x;
    x=x+s;
    fc=feval(f,x);
    itc=itc+1;
%
% Too many iterations?
%
    if(itc > maxit)
        disp(' maxit reached');
        if nargout == 2
            histmp(itc+1,:)=[itc, fc];
            hist=histmp(1:itc+1,:);
        end
        return
    end
%
    if nargout == 2
        histmp(itc+1,:)=[itc,fc];
    end
end
if nargout == 2
    hist=histmp(1:itc+1,:);
end