% Root finding using inverse interpolation
% a, b, & c are 3 initial guesses to the root.
% See p.234 of Heath, example 5.13
% K. Ming leung, 02/08/03

clear all;
fcn='x2m4sinx'; format long;
a=1; fa=feval(fcn,a);
b=2; fb=feval(fcn,b);
c=3; fc=feval(fcn,c);
itMax=20; it=1; tol=1e-8; h=10;
dsp=[];                    % for diaplay only
while it <= itMax & abs(h) > tol;
    u=fb/fc; v=fb/fa; w=fa/fc;
    p=v*(w*(u-w)*(c-b)-(1-u)*(b-a));
    q=(w-1)*(u-1)*(v-1);
    h = p/q;
    c=a; fc=fa;         % c replaced by the old a
    a=b; fa=fb;         % a replaced by the old b
    b=b+h; fb=feval(fcn,b);     % new b
    dsp=[dsp; [it b fb]]; it=it+1;
end;
display(dsp);