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;   returnend% set tolerances if necessaryif( nargin <= 3 ) tolf = 1.e-7; tolx = 1.e-7; maxit = 100; endif( nargin <= 4 ) tolx = 1.e-7; maxit = 100; endif( 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;endx = c; % check why the Regula falsi truncated and set iflagif( abs(b-a) >= tolx & abs(fc) >= tolf )    % the Regula Falsi truncated because the maximum number of iterations    % was reached    iflag = 1;    returnend