HELP! Bisection method problems

I've got a some problems with bisection algorithm. I paste here my code.

function [x, output] = fzer0v5_1_1(fun,xo,varargin)
narginchk(2, 5); % check if the function receives the right number of input parameters
nargoutchk(0,2); % check if the function receives the right number of output parameters

% PRIORITY CONTROL. IF BISECTION METHOD IS INAPPLICABLE OR FUN IS NOT A FUNCTION HANDLE => ERRORE IS GENERATE! 
% CHECK XO INTERVAL TOO. 
CheckXO(fun, xo); 

% OPTIONAL PARAMETER CONTROL
switch nargin
    case 5; [TOL, NMAX, g] = CheckG(varargin{1}, varargin{2}, varargin{3});
    case 4; [TOL, NMAX, g] = CheckNMAX(varargin{1}, varargin{2}, 'N');
    case 3; [TOL, NMAX, g] = CheckTOL(varargin{1}, 500, 'N');
    case 2; TOL = eps; NMAX = 500; g = 'N';
    otherwise
end
nit= 0;
XO=xo;
% FIND THE ZERO IF IT IS IN ONE OF THE TWO EXTREMES. IF WE HAVE 2 ZERO ON
% BOTH EXTREMES THEN, THE FUCNTION RETURN TO THE USER THE FIRST
a=find(~(arrayfun(fun,xo)),1); 
if ~isempty(a)
    x = xo(a);
else
% IF ZERO IS NOT FOUND ON EXTREMES, THEN THE MAIN LOOP OF ALGORITHM START
x = xo(1)+0.5*diff(xo);
TOLX = max(TOL * max(abs(xo)), realmin);
while(abs(diff(xo))>= TOLX && abs(fun(x))>=eps && nit<NMAX) 
    if ~(isequal(sign(fun(x)),sign(fun(xo(1))))) % THIS TYPE OF SIGN CHECK IS LESS ERROR PRONE (UNDERFLOW/OVERFLOW) THAN A COMMONPLACE MULTIPLICATION
        xo(2)= x;
    else
        xo(1)= x;
    end
    nit=nit+1;
    x = xo(1)+0.5*diff(xo);
    TOLX = max(TOL*max(abs(xo)), realmin);
end
end
% IF THE NMAX VALUE IS REACHED THEN THE USER WILL BE NOTIFIED WITH A WARNING
if(nit==NMAX) 
    warning("Limit reached... calculation stopped!");
end
if g == 'Y'
    plotFunction(fun, XO, x);
end
if (nargout ==2 )
    output.fx = fun(x);
    output.niter = nit;
end
end

% CHECK FUNCTION OF G PARAMTER
function [TOL, NMAX, g] = CheckG(tol, nmax, G)
if ~ismember(G, {'Y','N'})
    error("g charater must be Y or N. See doc fzer0");
end
[TOL, NMAX, g] = CheckNMAX(tol, nmax, G);
end

% CHECK FUNCTION OF NMAX PARAMTER
function [TOL, NMAX, g] = CheckNMAX(tol, nmax, G)
if (~isreal(nmax)||~isnumeric(nmax)||~(length(nmax)==1)||~isfinite(nmax)||nmax<=0 )
    error('NMAX must be a real number greather than 0.');
elseif ~(floor(nmax)==nmax)
    error("NMAX must be an integer value.");
end
[TOL,NMAX,g] = CheckTOL(tol,nmax, G);
end

% CHECK FUNCTION OF TOL PARAMETER
function [TOL, NMAX, g] = CheckTOL(tol,nmax, G)
if (~isreal(tol)||~isnumeric(tol)||~(length(tol)==1))
    error('TOL must be a real number.');
elseif ~isfinite(tol)
    error("TOL can't be + - Inf or NaN.");
elseif tol < eps || tol >=1 
    warning("tolerance value is wrong, see doc fzer0. Setted to default value.");
    tol = eps;
end
TOL = tol; NMAX = nmax; g = G;  
end

% CHECK FUNCTION OF FUN PARAMETER
function CheckFun(FUN,XO)
if ~ isa(FUN,'function_handle')
    error("First argoument must be a function handle of a in R function.");
elseif FUN(XO(1))*FUN(XO(2))>0  
    error("Method not applicable for this interval. See doc fzer0.")
end
end

% CHECK FUNCTION OF XO PARAMTER
% ATTENTION TO DEATILS => IF XO(1) IS GREATHER THAN XO(2) THIS FUNCTION IS HOWEVER APPLICABLE!
function CheckXO(FUN, XO) 
if (~isreal(XO)||~isnumeric(XO))
    error('xo must be composed by only real number.');
elseif ~all(isfinite(XO))
    error("xo elements can't be + - Inf or NaN.");
elseif (~isvector(XO)||~(length(XO)==2))
    error('xo must be a vector of 2 elements.');
end
CheckFun(FUN,XO);
end

I have problem when i must check the sign of the product of the function in "CheckFun". Multiplication is underflow/overflow prone and I have to write robust code. So I thought of using a solution like "if" in the while of the algorithm, but I have a problem in situation like:

>> sin(pi)

ans =

   1.2246e-16

>> sin(0)

ans =

     0

If I specified the interval [0,pi], my program is able to find the first zero, but that control in "CheckFun" returns an error if I test signs of both function value... How can I do to solve this problem?