function [b,fval,exitflag,output] = fzero(FunFcn,x,options,varargin) %#ok<STOUT>
%MATLAB Code Generation Library Function

%   Limitations:
%   1. The first argument must be a function handle.  Struct, inline
%   function, and string inputs for the first argument are not supported.
%   2. Up to 3 output arguments are supported. The 4th output argument, a
%   struct, is not supported.
%   3. Only two fields of an option struct are supported: TolX and
%   FunValCheck.  All other options in an option struct input are ignored.

%   Copyright 1984-2023 The MathWorks, Inc.
%#codegen

eml_invariant(nargin > 0, ...
    'MATLAB:minrhs');
eml_invariant(nargout < 4, ...
    'Coder:toolbox:fzero_2', ...
    'IfNotConst','Fail');
eml_invariant(~isstruct(FunFcn), ...
    'Coder:toolbox:OptimNoProblemStruct', upper(mfilename), ...
    'IfNotConst','Fail');

defaultopt = struct('Display',     'notify', ...
                    'MaxIter',     [], ...
                    'MaxFunEvals', [], ...
                    'TolX',        eps, ...
                    'TolFun',      [], ...
                    'FunValCheck', 'off', ...
                    'OutputFcn',   [], ...
                    'PlotFcns',    []);

% Undocumented syntax used by OPTIMSET to query default options structs
if nargin == 1 && eml_is_const(FunFcn) && strcmp(FunFcn, 'defaults')
    b = defaultopt;
    return;
end

eml_invariant(isa(FunFcn,'function_handle'), ...
    'Coder:toolbox:fzero_4', ...
    'IfNotConst','Fail');
eml_invariant(nargin >= 2, ...
    'MATLAB:optimfun:fzero:NotEnoughInputs', ...
    'IfNotConst','Fail');
eml_invariant(isa(x,'double'), ...
    'Coder:toolbox:fzero_6', ...
    'IfNotConst','Fail');
eml_invariant(eml_is_const(size(x)), ...
    'Coder:toolbox:OptimfunFixedSizeX', ...
    'IfNotConst','Fail');
eml_invariant(isscalar(x) || numel(x) == 2, ...
    'Coder:toolbox:fzero_8', ...
    'IfNotConst','Fail');


if nargin < 3 
    options = defaultopt;
else
    eml_invariant(isstruct(options) || (eml_is_const(isempty(options)) && isempty(options)), ...
                  'Coder:toolbox:OptimfunOptArg', ...
                  'IfNotConst','Fail');
end

tol = double(optimget(options,'TolX',defaultopt,'fast'));
funValCheck = strcmp(optimget(options,'FunValCheck',defaultopt,'fast'), 'on');

% Initialization
exitflag = 1;
% Add a wrapper function to check for Inf/NaN/complex values
% Interval input
if isscalar(x)
    eml_invariant(isfinite(x), ...
        'MATLAB:optimfun:fzero:Arg2NotFinite');
    fx = FunFcn(x,varargin{:});
    eml_invariant(isreal(fx), ...
        'Coder:toolbox:fzero_11', ...
        'IfNotConst','Fail');
    if fx == 0
        b = x;
        fval = fx;
        return
    elseif ~isfinite(fx)
        coder.internal.error('MATLAB:optimfun:fzero:ValueAtInitGuessComplexOrNotFinite');
    end
    if x ~= 0
        dx = eml_rdivide(x,50);
    else
        dx = 0.02;
    end
    % Find change of sign.
    twosqrt = sqrt(2);
    a = x; fa = fx; b = x; fb = fx;
    while (fa > 0) == (fb > 0)
        dx = twosqrt*dx;
        a = x - dx;  fa = FunFcn(a,varargin{:});
        if funValCheck
            checkval(fa);
        end
        if ~isfinite(fa)
            exitflag = -3;
            b = coder.internal.nan; fval = coder.internal.nan;
            return
        elseif ~isfinite(a)
            exitflag = -6;
            b = coder.internal.nan; fval = coder.internal.nan;
            return
        end
        if (fa > 0) ~= (fb > 0) % check for different sign
            break
        end
        b = x + dx;  fb = FunFcn(b,varargin{:});
        if funValCheck
            checkval(fa);
        end
        if ~isfinite(fb)
            exitflag = -3;
            b = coder.internal.nan; fval = coder.internal.nan;
            return
        elseif ~isfinite(b)
            exitflag = -6;
            b = coder.internal.nan; fval = coder.internal.nan;
            return
        end
    end % while
    savefa = fa; savefb = fb;
else
    eml_invariant(isfinite(x(1)) && isfinite(x(2)), ...
        'MATLAB:optimfun:fzero:Arg2NotFinite');
    a = x(1);
    b = x(2);
    fa = FunFcn(a,varargin{:});
    eml_invariant(isreal(fa), ...
        'Coder:toolbox:fzero_14', ...
        'IfNotConst','Fail');
    fb = FunFcn(b,varargin{:});
    if ~(isfinite(fa) && isfinite(fb))
        coder.internal.error('MATLAB:optimfun:fzero:ValuesAtEndPtsComplexOrNotFinite');
    end
    savefa = fa; savefb = fb;
    if fa == 0
        b = a;
        fval = fa;
        return
    elseif fb == 0
        % b = b;
        fval = fb;
        return
    elseif (fa > 0) == (fb > 0)
        coder.internal.error('MATLAB:optimfun:fzero:ValuesAtEndPtsSameSign');
    end
    % Starting guess scalar input
end % if isscalar(x)
fc = fb;
% The following initializations are required by the compiler.  All will be
% overwritten at the top of the loop.
c = b;
e = 0;
d = 0;
% Main loop, exit from middle of the loop
while fb ~= 0 && a ~= b
    % Ensure that b is the best result so far, a is the previous
    % value of b, and c is on the opposite side of the zero from b.
    if (fb > 0) == (fc > 0) % Always evaluates to true on first iteration.
        c = a;  fc = fa;
        d = b - a;  e = d;
    end
    if abs(fc) < abs(fb)
        a = b;    b = c;    c = a;
        fa = fb;  fb = fc;  fc = fa;
    end
    % Convergence test and possible exit
    m = 0.5*(c - b);
    toler = 2.0*tol*eml_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 = eml_rdivide(fb,fa);
        if (a == c)
            % Linear interpolation
            p = 2.0*m*s;
            q = 1.0 - s;
        else
            % Inverse quadratic interpolation
            q = eml_rdivide(fa,fc);
            r = eml_rdivide(fb,fc);
            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 = eml_rdivide(p,q);
        else
            d = m;  e = m;
        end
    end % Interpolation
    % Next point
    a = b;
    fa = fb;
    if abs(d) > toler
        b = b + d;
    elseif b > c
        b = b - toler;
    else
        b = b + toler;
    end
    fb = FunFcn(b,varargin{:});
    if funValCheck
        checkval(fb);
    end
end % Main loop
fval = fb; % b is the best value
if abs(fval) <= eml_max(abs(savefa),abs(savefb))
else
    exitflag = -5;
end

%------------------------------------------------------------------

function checkval(x)
% Note: we do not check for Inf as FZERO handles it naturally.  ???
if isnan(x)
    coder.internal.error('Coder:MATLAB:fzero_checkfun_NaNFval',coder.internal.flt2str(x));
end

%--------------------------------------------------------------------------

% LocalWords:  Evals Pts checkfun NFval
