bisect := proc(  f :: algebraic, 
                 x :: name,
                 a :: numeric, 
                 b :: numeric, 
               tol :: numeric)
   local aa, bb, mp;
 
   aa := evalf(a); bb := evalf(b);  # transfer a, b

   if aa > bb then
      print("a must be smaller than b");
      return;
   end if;

   if evalf(subs(x=aa,f)*subs(x=bb,f) ) > 0 then
      print("The function must have different sign at a and b" );
      return;
   end if;
   
   while bb - aa > tol do

      mp := evalf((aa+bb)*0.5);  # the mid-point of [aa,bb]
      if evalf(subs(x=aa,f)*subs(x=mp,f)) <= 0 then
         bb := mp;    # update the right end
      else
         aa := mp;    # update the left end
      end if;

      print( evalf((aa+bb)/2), abs(bb-aa));

   end do;

   return evalf((aa+bb)/2) ;  # output the approx. zero

end proc;