function [a, b, x, ithist, iflag] = regulafalsi( f, a, b, tolf, tolx, maxit )
%
%  function [a, b, x, ithist, iflag] = regulafalsi( f, a, b, tolf, tolx, maxit )
%
%  use Regula Falsi to compute an approximate root of f.
%
%  Input parameters:
%    f       name of a matlab function that evaluates 
%    a, b    real numbers satisfying f(a)*f(b) < 0
%    tolf    stopping tolerance (optional. Default tolf = 1.e-7)
%            Regula Falsi stops if |f(x)| < tolf
%    tolx    stopping tolerance (optional. Default tolx = 1.e-7)
%            Regula Falsi stops if  b-a < tolx
%    maxit   maximum number of iterations (optional. Default maxit = 100)
%
%
%  Output parameters:
%    a, b    real numbers satisfying f(a)*f(b) < 0; a root of f is
%            located between a and b
%    x       approximation of the solution. 
%    ithist  array with the iteration history
%            The i-th row of ithist contains  [it, a, b, c, fc]
%    ifag    return flag
%            iflag = -1  error in input data, f(a)*f(b) > 0
%            iflag =  0  |b-a| < tolx and |x-x*| < 0.5*tolx, where x*
%                        is a root of f
%            iflag =  1  iteration terminated because maximum number of 
%                        iterations was reached. |b-a| >= tolx
%
%
%  Matthias Heinkenschloss
%  Department of Computational and Applied Mathematics
%  Rice University
%  Jan 17, 2002
%
%

fa = feval(f, a);
fb = feval(f, b);
if( fa*fb > 0 );
   iflag = -1;
   return
end


% set tolerances if necessary
if( nargin <= 3 ) tolf = 1.e-7; tolx = 1.e-7; maxit = 100; end
if( nargin <= 4 ) tolx = 1.e-7; maxit = 100; end
if( nargin <= 5 ) maxit = 100; end
   
fc = 2*tolf;


it    = 0;
iflag = 0;
while( it < maxit & abs(b-a) >= tolx & abs(fc) >= tolf )
   c  = a - fa * (b-a)/(fb-fa);
   fc = feval(f, c);
   
   ithist(it+1,:) = [it, a, b, c, fc];
   
   if( fa*fc < 0 );  
       b = c; fb = fc;
   else
       a = c; fa = fc;
   end
   it = it+1;
end

x = c; 

% check why the Regula falsi truncated and set iflag
if( abs(b-a) >= tolx & abs(fc) >= tolf )
    % the Regula Falsi truncated because the maximum number of iterations
    % was reached
    iflag = 1;
    return
end