function [c] = bisect(f, x0, x1, tolerance, maxiter, printflag)
% bisect find root of f(x) = 0 using bisection method
% f = function handle of function whose root is desired
% x0 = left end of initial interval to be bisected
% x1 = right hand end of initial interval
% tolerance = desired acuracy of root
% maxiter = maximum bisections to try
% printflag = true, show intermediate results
% printflaf = false, show only the final answer
% c = best approximation to root is returned

   a = x0;
   fa = f(a);
   b = x1;
   fb = f(b);
   
   if sign(fa) == sign(fb)
      error('function has same sign at end points');
   end
   
   if printflag
      fprintf('%3s%19s%19s%19s%19s\n', 'n','a','c','b','err');
   end

   n = 1;
   while n <= maxiter
   
      % bisect interval [a,b] and calculate function
      % value at mid-point
   
      err = 0.5*(b-a);
      c = a + err;
      fc = f(c);
   
      if fc == 0 || err < tolerance
         break;
      end
   
      if printflag
         fprintf('%3d%19.10e%19.10e%19.10e%19.10e\n', n,a,c,b,err);
      end
     
      % choose next subinterval [a,b] using bisection
      
      if sign(fa) ~= sign(fc)
         % root is in left half [a,c] of [a,b]
         b = c; 
         fb = fc;
      elseif sign(fc) ~= sign(fb)
         % root is in right half [c,b] of [a,b]
         a = c;
         fa = fc;
      end
   
      n = n + 1;   
   end
   
   if n > maxiter
      fprintf('No convergence for max iterations %d\n', n-1);
   end
end