function [x,num_its] = bisect_td(f,a,b,tol)
% bisection method
% check error conditions
if (f(a)*f(b) >= 0.0)
  error('f(a)*f(b)>=0.0');
end
if (a>=b)
  error('a >=b');
end
ak=a;
bk=b;
num_its=0;

while(bk - ak > 2*tol)
    c=0.5*(bk+ak);
  if (f(c)*f(bk)>=0.0)
% x is in [ak,c]       
    bk = c;
  else  
    ak = c;
  end % if
  num_its = num_its+1;
end  % while  
x = 0.5*(ak+bk);  % answer

end

