function b = fzero2(xs,fx,c,erc)          

% function x=fzero2(xs,fx,c,erc)  find zero of one variable function like ssec2
%
% This program is a modified version of fzero for which the
% arguments are entered in much the same way as those for ssec2.
%
% xs gives an interval to seek a solution to f(x) = 0 in.
% fx is a character vector that tells how to find f(x).
% c is a vector of parameters in f(x)
% erc gives the difference between successive values of x allowed
%     when convergence is deemed to be achieved.
% If erc is not given, then the program uses:
%     erc=eps (eps is Matlab epsilon)
% If only two arguments are given, c is set to null.
% 
% Example 1: >> cs=[1 -2 -4.25 7.5];
%            >> disp(fzero2([2 3],'polyval(c,x)',cs,.0001))
% Example 2: >> disp(fzero2([3 4],'sin(x)'))
%
% ACP, DAJ

% Initialization
if nargin < 4, erc = eps; end
if nargin == 3, par=[]; end
a=xs(1);
x=a;
fa=eval(fx);
b=xs(2);
x=b;
fb=eval(fx);

fc1 = fb;
% Main loop, exit from middle of the loop
while fb ~= 0
   % Insure that b is the best result so far, a is the previous
   % value of b, and c1 is on the opposite of the zero from b.
   if (fb > 0) == (fc1 > 0)
      c1 = a;  fc1 = fa;
      d = b - a;  e = d;
   end
   if abs(fc1) < abs(fb)
      a = b;    b = c1;    c1 = a;
      fa = fb;  fb = fc1;  fc1 = fa;
   end

   % Convergence test and possible exit
   m = 0.5*(c1 - b);
   toler = 2.0*erc*max(abs(b),1.0);
   if (abs(m) <= toler) + (fb == 0.0), break, end

   % Choose bisection or interpolation
   if (abs(e) < toler) + (abs(fa) <= abs(fb))
   % Bisection
      d = m;  e = m;
   else
   % Interpolation
      s = fb/fa;
      if (a == c1)
      % Linear interpolation
         p = 2.0*m*s;
         q = 1.0 - s;
      else
      % Inverse quadratic interpolation
         q = fa/fc1;
         r = fb/fc1;
         p = s*(2.0*m*q*(q - r) - (b - a)*(r - 1.0));
         q = (q - 1.0)*(r - 1.0)*(s - 1.0);
      end;
      if p > 0, q = -q; else p = -p; end;
      % Is interpolated point acceptable
      if (2.0*p < 3.0*m*q - abs(toler*q)) * (p < abs(0.5*e*q))
         e = d;  d = p/q;
      else
         d = m;  e = m;
      end;
   end % Interpolation
   % Next point
   a = b;
   fa = fb;
   if abs(d) > toler, b = b + d;
   else if b > c1, b = b - toler;
        else b = b + toler;
        end
   end
   x=b;
   fb=eval(fx);

end % Main loop
