function [x, rate] = bisect( f, a, b, tol, maxIter, verbose = false )
  % ----------------------------------------------------------------------
  % Jonathan Senning <jonathan.senning@gordon.edu>
  % Gordon College
  %
  % Usage: x =  bisect(@f, a, b, tol, maxIter[, verbose])
  %        or
  %        [x,rate] = bisect(@f, a, b, tol, maxIter[, verbose])
  %
  % This function uses bisection to solve f(x) = 0.
  %
  % This function requires that the function f(x) be already defined.
  % ----------------------------------------------------------------------

  if ( nargin ~= 5 && nargin ~= 6 )
    fprintf( '\aUsage: ' );
    fprintf( '[x,rate] = bisect(@f, a, b, tol, maxIter[, verbose])\n' );
    return
  end

  x  = a;
  fa = f( a );
  fb = f( b );

  % Make sure that (according to the Intermediate Value Theorem) the
  % specified interval does contain a root.

  if ( sign( fa ) == sign( fb ) )
    fprintf( 'Interval [%f,%f] may not contain a root\n', a, b );
  end

  % Initialize for the search.  Note that we try and minimize the number
  % of evaluations of f(x).  It is initially evaluated twice before the
  % iteration begins and then once during each iteration. These extra
  % points are needed to start off the error calculation which compares
  % estimates of solutions just computed with previous estimates.

  x0 = x + 4.0 * tol;
  x1 = x + 8.0 * tol;

  k = 0;

  if verbose
    fprintf( '     k         a                  b          ' );
    fprintf( ' Estimate of root  Max. error\n' );
    fprintf( '    -- ------------------ ------------------ ' );
    fprintf( '------------------ ----------\n' );
  end

  while ( k <= maxIter )

    x2 = x1;
    x1 = x0;
    x0 = x;
    x = ( a + b ) / 2.0;
    fx = f( x );

    max_error = abs( b - a ) / 2.0;

    if verbose
      fprintf( '    %2d %18.11e %18.11e %18.11e %10.4e\n', ...
	       k, a, b, x, max_error );
    end

    if ( max_error < tol )
        break;
    end

    k = k + 1;

    if ( sign( fa ) == sign( fx ) )
      a = x;
      fa = fx;
    else
      b = x;
    end

  end

  if ( k > maxIter )
    fprintf( 'Error: exceeded %d iterations\n', maxIter );
  end

  rate = log( abs((x - x0) / (x0 - x1)) ) / log( abs((x0 - x1) / (x1 - x2)) );

end
