% Root finding using linear fractional interpolation
% a, b, & c are 3 initial guesses to the root.
% See p.235 of Heath, example 5.14
% K. Ming Leung, 03/13/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=[];
while it <= itMax & abs(h) > tol;
    ac=a-c; bc=b-c;
    h = fc*ac*bc*(fa-fb)/(fa*ac*(fc-fb)-fb*bc*(fc-fa));
    a=b; fa=fb;     % a replaced by old b
    b=c; fb=fc;     % b replaced by old c
    c=c+h; fc=feval(fcn,c);     % new c
    dsp=[dsp; [it c fc]]; it=it+1;
end;
display(dsp);