#
# Bisection method for solving f(x)=0
#
#        Input:   f -- function. 
#                         (e.g. [> f := x -> x^2-3; )
#               a,b -- the interval where f(a)*f(b)<0
#               tol -- the accuracy tolerance
#     
bisec := proc(    f :: operator,
                  a :: numeric, 
                  b :: numeric,
                tol :: numeric
                                 )
   local aa, bb, mp;

   aa, bb := a, b;  # pass the first interval

   if evalf(f(aa)*f(bb)) > 0 then
      print(`There may not be a solution`);
      RETURN();
   end if;

   printf(`    a              b    `);
   printf(`          mid-point         error \n`);
   while  abs(bb-aa) > tol  do

      mp := 0.5*(aa+bb); # update the mp

      printf(`%14.9f,%14.9f,%14.9f,     %8.2e \n`,
              aa,    bb,     mp,    bb-aa);

      if evalf(f(aa)*f(mp)) < 0 then
         
         bb := mp
  
      else
  
         aa := mp

      end if;

   end do;

   print( mp );

end proc;