% Bisection method for the function f based on the interval [a,b] containing
% the true root.
% Returns all the iterates in addition to the number of iterations.
function [x,iter] = bisection_method(a,b,f,ep)

MAX_ITER = 1000;
fva = f(a);
k = 0;
x = zeros(MAX_ITER,1);
while ( b-a > 2*ep & k < MAX_ITER )
   % Bisect the interval
   x(k+1) = (b+a)/2;
   % Get the new value of the function at the bisection points
   fvb = f(x(k+1));
   % Check if a root was found.
   if ( fvb == 0 )
      return
   end
   % Determine the new interval.
   if (fva*fvb < 0 )
      b = x(k+1);
   else
      a = x(k+1);
      fa = fvb;
   end
   % Increment the counter.
   k = k+1;
end
% Bisect one more time to ensure the root is within epsilon
k = k+1;
x(k) = (b+a)/2;
% Trim off the excess zeros.
x = x(1:k);

% If the iteration proceeded beyond MAX_ITER then display a warning message.
if k > MAX_ITER
   warning('Maximum number of iterations exceeded');
end
iter = k;