clear all; close all;
xp=0:0.1:5;, yp=f1(xp); plot(xp,yp)
axis([0.5 3.2 -5 22.5]); hold
yp0=0*xp; plot(xp,yp0,':')

x=3; tol=1e-12; plot(x,f1(x),'*')
x2=1;

% Tangent line: y-yp = f'(x) (x-xp)
k = length(xp);
for i=1:k
  yt(i) = f1(x) + (f1(x)-f1(x2))*(xp(i)-x)/(x-x2);
end
plot(xp,yt,':')

n = 0; nmax= 500; 
while abs(f1(x)) > tol &&  n< nmax
    n = n+1;
    % x=x-f1(x)/df1(x)
     x3=x-f1(x)*(x-x2)/(f1(x)-f1(x2))
    plot(x3,f1(x3),'o')
  
    for i=1:k
      yt(i) = f1(x) + (f1(x)-f1(x3))*(xp(i)-x)/(x-x3);
    end
    x2=x; x=x3;
    if n==1 
        plot(xp,yt,':')
    end
end