#
# program of bisection method solving
#           f(x) = 0
#
#  Input      f --- the function
#                     (e.g.   f:=x->x^2-3 )
#           a,b --- the interval [a,b]
#           tol --- the error tolerance
#
bisect:=proc(f,a,b,tol)
   local aa, bb, mp;

   if evalf( f(a)*f(b) ) > 0  then
      print(`invalid interval`);
   else
      aa:=a; bb:=b;
      while abs(bb-aa)>tol do
         mp:=0.5*(aa+bb);
         if evalf( f(mp)*f(aa) ) > 0 then
            aa:=mp;
         else
            bb:=mp;
         fi;
      od;
   fi;

   print(`The approximate solution`);
   0.5*(aa+bb);

end;