function root = bisect ( f, a, b )

fa = feval ( f, a ); 
fb = feval ( f, b );

if ( sign ( fa ) == sign ( fb ) )
  'BISECT - Fatal error!'
  '  [A,B] is not a change-of-sign interval.'
  root = 0;
  return
end

format long

while ( abs ( b - a ) > 0.000001 )

  c = ( a + b ) / 2;
  fc = feval ( f, c );

  [  a,  c,  b ;
    fa, fc, fb ]

  if ( fc == 0 ) 
    root = c;
    return
  elseif ( sign ( fc ) == sign ( fa ) )
    a = c;
    fa = fc;
  elseif ( sign ( fc ) == sign ( fb ) )
    b = c;
    fb = fc;
  end

end

root = 0.5 * ( a + b );