function [ x, a, b, step_num ] = false_position ( xtol, ftol, a, b, step_max, f )

%% false_position solves f(x)=0 by the method of false position.
%
% xtol, ftol, tolerance for interval and function size.
% a, b, the endpoints of a change of sign interval.
% step_max, maximum number of steps to take.
%
% x, output, the estimated solution.
% a, b, output, the updated endpoints of the change of sign interval.
% step_num, the number of bisection steps.
%
  if ( sign ( f(a) ) * sign ( f(b) ) ~= -1 )
    error ( 'F(A) and F(B) are not of opposite signs!' );
  end

  step_num = 0;

  while ( xtol < b - a )

    if ( step_max <= step_num )
      fprintf ( 'Maximum number of steps exceeded!\n' );
      return
    end

    step_num = step_num + 1;
%
%  Use the secant approximation.
%
    x = a - f(a) * ( b - a ) / ( f(b) - f(a) );
%
%  Is the function value small enough?
%
    if ( f(x) < ftol )
      return
    end
%
%  Choose new interval [A,B=X] or [A=X,B].
%
    if ( sign ( f(x) ) == sign ( f(a) ) )
      a = x;
    else
      b = x;
    end

  end

  return
end
