% based on Michael Heath, Scientific Computing, Ch.5

tic

itmax=100;

k=0;

fa=f(a);
b=a-fa/Df(a); k=k+2;

fb=f(b);
c = b + (b-a)/(fa/fb-1); k=k+1;

k
[a b c]



while abs(c-b)>tol*max(abs(c),1.0)
     if k+1>itmax
       break
     end 
     fc=f(c);
     u=fb/fc;
     v=fb/fa;
     w=fa/fc;
     p=w*(u-w)*(c-b)-(1-u)*(b-a);
     p=v*p;     
     q=(u-1)*(v-1)*(w-1);
     x=b+p/q;
     k=k+1
     a=b;
     fa=fb;
     b=c;
     fb=fc;
     c=x;
     [c abs(c-b)]
end  

toc