clear

f=inline('exp(x)+1');
x=-2:.01:2;
figure;plot(x,f(x));hold on;plot([x(1) x(end)],[0 0])

x=[-1 0 1];count=0;
while abs(f(x(3)))>1e-8 && count<100
    a=((x(2)-x(3))*(f(x(1))-f(x(3)))-(x(1)-x(3))*(f(x(2))-f(x(3))))/...
        ((x(1)-x(2))*(x(2)-x(3))*(x(1)-x(3)));
    b=((x(1)-x(3))^2*(f(x(2))-f(x(3)))-(x(2)-x(3))^2*(f(x(1))-f(x(3))))/...
        ((x(1)-x(2))*(x(2)-x(3))*(x(1)-x(3)));
    c=f(x(3));
    
    roots=x(3)+(-b+[1 -1]*sqrt(b^2-4*a*c))/(2*a);
    [minval,ind]=min(abs(roots-x(3)));
    newx=roots(ind)
    x=[x(2) x(3) newx];
    
    count=count+1;
end
x(3)