function [x,ierr] = bisection(f, a, b, tol, tolf)
%
% use bisection to calculate a root of f. 
%
% calling sequence: [x,ierr] = bisection(@f, a, b, tol)
%
% f = name of the matlab function to evaluate function f
% a,b = initial endpoints of the interval. Must have a<b, with
%       f(a) and f(b) different signs (otherwise return ierr=1)
% tol = terminate when length of interval falls below this value
% tolf = terminate if |f(x)| falls below this value, where x is some
%        point evaluated during the algorithm
% 
% return x    = final estimate of the root.
%        ierr =  error code (=0 if an approximate root is found)
%