ode15s unable to meet integration

Hello,

I am modeling a continously stirred reactor and some chemical reactions and trying to do a cyclic process, i.e. feed the reactive from 0 to a certain time t1 (I call this the reactive phase) then feed an inert from time t1 to a time t2 (I call this the inert phase) and so on. I’ve written the differential equations in a function called fun. I can’t provide the whole code because it involves a lot of matlab functions but the general idea is the following:

function [dydt] = fun(t,y)
N=250;
ygas=y(1:N);
yother=y(N+1:end);
 
%% ~~~~~~~~~~~~~~~ gas components ~~~~~~~~~~~~~~~
dygasdt=fgas(ygas);
 
%% ~~~~~~~~~~~~~~~ other ~~~~~~~~~~~~~~~
dyotherdt=fother(yother,ygas);
 
%%
dydt=[dygasdt;dyotherdt];
end

I solve with ode15s:

% the time vectors and initial conditions are defined (tps_reac and tps_inert are the vectors indicating the final time of each phase)
for i=1:ncycles
    %% reactive is fed 
    [t,y]=ode15s(@(t,y)fun(t,y),0:1:tps_reac(i),yo);
    yo=y(end,:);
    %% inert is fed 
    [t1,y1]=ode15s(@(t,y)fun(t,y),0:1:tps_inert(i),yo);
    yo=y1(end,:);   
%then comes the part of the code where the results are stored in respective time and solution matrices
end

The initial conditions for the first reactive phase are set as a molar fraction of 1 for the reactive (and 0 for the other components). The initial condition for the first inert phase is set as the solution from the reactive phase at tps_reac(1), the initial condition for the second reactive phase is set as the solution from the first inert phase at tps_inert(1) and so on.

My problem is the following, during the inert phase what is fed into the reactor is a pure inert species and so the molar fraction of the inert species quickly reaches a value close to 1 and the molar fraction of the other gas species quickly reach a value close to 0, sometimes they even become negative (it’s an open reactor so the species are pushed out of the reactor by the inert gas). So after all the mixture that was left over from the reactive phase is pushed out, nothing much happens in the gas phase. However, the “other” phase is still calculated. I get the following type of error:

Warning: Failure at t=1.927609e+01.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.684342e-14) at time t.

I am guessing it’s because the orders of magnitude calculated for the gas phase become too small as the integration goes on. Is there a way to tell the code to not calculate dygasdt or to set ygas to 0 for values of ygas that are lower than 1e-15 for example? Would this help the integration?