The Bisection Method

function [a,b,fail] = bisect(f,a,b,eps,N)
%
%  Applies the bisection method to the function f on [a,b],
%  f is assumed to be continuous and f(a)*f(b) < 0.
%  Stops when |f(a)| < eps or |f(b)| < eps or |b-a| < eps.
%  Failure is reported (fail = 1) when the accuracy requirement
%  is not satisfied in N steps; otherwise fail = 0 on return.
%
%  Example : 
%  >> [a,b,fail] = bisect('cos',pi/4,2*pi/3,0.0001,100)
%
fprintf('running the bisection method...\n');
for i = 1:N 
   m = (a+b)/2;
   fm = feval(f,m);
   if fm*feval(f,a) < 0
      b = m;
   else
      a = m;
   end;
   fprintf('  a = %f, b = %f, m = %f, f(m) = %f\n', a,b,m,fm);
   if (abs(fm) < eps) | ((b-a) < eps)
      fail = 0;
      fprintf('succeeded after %d steps\n', i);
      return;
   end
end
fprintf('failed requirements after %d steps\n', N);
fail = 1;

The output of this algorithm on cos(x) = 0 with [a,b] = [pi/4,2*pi/3] :

>> [a,b,fail] = bisect('cos',pi/4,2*pi/3,0.0001,100);
running the bisection method...
  a = 1.439897, b = 2.094395, m = 1.439897, f(m) = 0.130526
  a = 1.439897, b = 1.767146, m = 1.767146, f(m) = -0.195090
  a = 1.439897, b = 1.603521, m = 1.603521, f(m) = -0.032719
  a = 1.521709, b = 1.603521, m = 1.521709, f(m) = 0.049068
  a = 1.562615, b = 1.603521, m = 1.562615, f(m) = 0.008181
  a = 1.562615, b = 1.583068, m = 1.583068, f(m) = -0.012272
  a = 1.562615, b = 1.572842, m = 1.572842, f(m) = -0.002045
  a = 1.567728, b = 1.572842, m = 1.567728, f(m) = 0.003068
  a = 1.570285, b = 1.572842, m = 1.570285, f(m) = 0.000511
  a = 1.570285, b = 1.571563, m = 1.571563, f(m) = -0.000767
  a = 1.570285, b = 1.570924, m = 1.570924, f(m) = -0.000128
  a = 1.570605, b = 1.570924, m = 1.570605, f(m) = 0.000192
  a = 1.570764, b = 1.570924, m = 1.570764, f(m) = 0.000032
succeeded after 13 steps