function [x, hist] = chordsol(x, f, tola, tolr)
%CHORDSOL  Chord method code for scalar equations.
%
% This code terminates on small relative-absolute errors
%
% function [X, HIST] = CHORDSOL(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 this out by calling chordsol as
%                Z = CHORDSOL(X, F, TOLA, TOLR) 
%                the storage for hist is not allocated and
%                the iteration history is not stored.
%
% Uses: fordiff.m, Tim's numerical derivative
%
%
maxit=20;    % maximum iterations
alpha=1.d-4; 
small=1.d-8;
if nargout == 2
    histmp=zeros(maxit,2);
end
%
% Initialize
%

itc=0;
fc=feval(f,x);
fp=fordiff(f,x,fc);
tol=tolr*abs(fc)+tola;
%
%    Store iteration history?
%
if nargout == 2
    histmp(itc+1,:)=[itc, abs(fc)];
end
%
%    Main chord 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)];
            hist=histmp(1:itc,:);
        end
    return
    end
%
%    Check for small relative derivatives.
%
    if(abs(fp) <= small*abs(fc))
        disp(' small derivative error');
    return;
    end
    s=-fc/fp;
    xt=x+s;
    ft=feval(f,xt);
    x=xt;
    itc=itc+1;
%
%    Too many iterations?
%
        if(itc > maxit)
            disp(' maxit reached');
            if nargout == 2
                histmp(itc+1,:)=[itc,abs(fc)];
                hist=histmp(1:itc,:);
            end
            return
        end
%
        fc=feval(f,x);
    if nargout == 2
        histmp(itc+1,:)=[itc,abs(fc)];
        end
    end
        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;