function c = Bisection(f, a, b, n_max, epsilon)
% Bisection root finding function 
% Where f is a function, and a and b are endpoints of the inverval we are searching.
% n_max is the maximum cycles we will attempt.
% epsilon is the tolerance for determining if we have found a root.

fa = f(a); 
fb = f(b); 
fprintf('interval [%6.1f, %6.1f]', a,b) 
fprintf('\nf(a) = %8.1f, f(b) = %8.1f\n', fa, fb) 

if sign(fa)==sign(fb) 
   disp('No root in the interval') 
   c = 0; 
else 
   error = b-a; 
   for n=0:n_max 
      error = error/2; 
      c = a + error; 
      fc = f(c); 
      fprintf('\n%6d, c = %10.8f fc = %10.8f error = %10.8f ', n, c, fc, error) 

      if error < epsilon 
           fprintf('\nconvergence\n') 
           break; % root was found 
      else 
           if sign(fa) ~= sign(fc) 
              b = c; 
              fb = fc; 
           else 
              a = c; 
              fa = fc; 
           end 
       end 
   end 
end