clear all;
format long

xm1=1; x0=2; func='x^2-2'; nmax=50; toll=10^-10; 

x=xm1; f=eval(func); fx=[f]; x=x0; f=eval(func); fx=[fx, f]; 
nit=0;     
xvect=[xm1,x0]; xdif=[]; f=toll+1; kprime=1;
while (nit < nmax) & (abs(f) > toll),
nit=nit+1; dim=length(xvect);
x=xvect(dim); fxk=eval(func); xk=x; i=dim;
while (i >= kprime), i=i-1; x=xvect(i); fxkpr=eval(func);
if ((fxkpr*fxk) < 0), xkpr=x; kprime=i; break; end;
end;
x=xk-fxk*(xk-xkpr)/(fxk-fxkpr); xvect=[xvect, x]; f=eval(func);
fx=[fx, f]; err=abs(x-xkpr); xdif=[xdif, err];
end;
xvect
format short


