% This is a simple implementation of Algorithm2.1, approx. 35 lines

% Declare the relevant parameters

tol=0.0001;
a=1.0;
b=2.0;
nmax=100;
p=2.0;
false = 1;
true = 0;
finish = false;
iteration = 0;

while (finish == false & iteration <= nmax );

   iteration = iteration + 1;
   pold=p;
%   p=(a+b)/2;
    p=a + (b-a)/2;

   fa=f(a);
   fb=f(b);
   fp=f(p);
   if (abs(p-pold)/abs(p) < tol) | (fp==0.0) ;

      % finished !

      p_final = p;
      finish = true;

   else;

      % not converged ! calculate new a and b

      if (fa*fp > 0);

         %  root is between p and b   
         a=p;

      else;

         %  root is between a and p
         b=p;

      end
    end

end

if (iteration < nmax);

  fprintf(1,'result after iteration number %5d ',iteration);
  fprintf(1,'is %15.7e \n',p_final);

else

  fprintf(1,'no convergence after iteration number %3d ',nmax);

end