How to solve a system of ODEs with a parameter discontinuity from min() function

Hello, I have a serious problem. I have a model that is positive for all values, but ODE45, ode15s etc, solve it and give results that are not correct. I think think that my time derivatives are discontinuous functions because of the min() function in the system. For certain parameter values, the results show that some of the state variables, could be negative. For instance, C, in the code below, takes negative at steady state, but mathematical if C happens to be zero, dc/dt will be increasing and thus we should always get positive values. But that is not teh result as you could see by running the code below. Could you please suggest how to resolve this.

Thank you very much for the anticipated assistance.

function grazingsteady mu_B=2;%2; epsilon=0; K_f=3e-6; K_g=15.9e-6;%C:N here is 5.3:1 normally in the pools it is 4-10:1 and in compost 20:30:1 mu_G=0.2; K_h=8e-3; %K_h=8e-3; alpha_mo=0.5; w=0.1*alpha_mo*mu_G;%0.01 r=0.5; D=0.001; %(0.25-0.7 e-3 per hr).01D not below 0.0001, for D=0.001, only G goes to 0 theta_G=1/10; theta_B=0.2;%1/5; C_in=159; N_in=30; time=[0,100]; %IC=[0.01,0.01,5.3*0.01,0.01]; IC=[0.5e-3,0.1e-3,2e-3,10e-6]; % GrowthB=mu_B*(IC(3))/(K_g+IC(3)) % cond2=theta_B-theta_G*alpha_mo % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % lambdaC=-(D*K_g)/(D-mu_B)% ensure that D<mu_B % Bac=r*(C_in-lambdaC) % Nitro=N_in-theta_B*r*(C_in-lambdaC) % condlHs=N_in/(C_in-lambdaC) % condRHS=theta_B*r %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % lambdaN=-(D*K_g)/(D-mu_B) % Bac1=(C_in-lambdaC)/theta_B % Carbon=C_in-(N_in-lambdaN)/(r*theta_B)% ensure that lambda_N<N_in and (N_in-labdaN)<rtheta_B % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % delta1=-((K_h*(D+w))/(-alpha_mo*mu_G+D+w))%ensure that D+w<alphamu_G % C1=(1/2*D*r)*sqrt(delta1^2*mu_B^2-2*delta1*C_in*D*r*mu_B+2*delta1*D*K_g*r*mu_B... % +C_in^2*D^2*r^2+2*C_in*D^2*K_g*r^2+D^2*K_g^2*r^2-delta1*mu_B+C_in*D*r-D*K_g*r) % G1=(delta1*(-D+mu_B*(C1/(K_g+C1)))/(mu_G*(delta1/(K_h+delta1)))) % N1=1/D*(D*N_in+theta_G*w*G1+(theta_B-theta_G*alpha_mo)*mu_G*(delta1/(delta1+K_g)*G1-theta_B*mu_B*delta1*(C1/(C1+K_g))))

% Di=[0:0.00025:0.001]; % De=zeros(1,length(Di)); % for i=1:length(Di) % D=Di(i); [t,y]=ode15s(@system,time,IC); % %[T,Z]=ode45(@system1,time,IC1); %Z(:,4) % De(i)=mean(y(:,5)); % set(0,'DefaultAxesFontSize',18,'DefaultTextFontSize',18); % plot(Di,De,'linewidth',2); % xlabel('Dilution rate (per Day)'); % ylabel('Mean degradation rate'); %plot(t,y(:,5),'r','linewidth',2); % end % % figure %set(0,'DefaultAxesFontSize',18,'DefaultTextFontSize',18); % plot(Di,De,'linewidth',2); % xlabel('Time(Days)'); % ylabel('Mean degradation rate'); %mean(y(:,5)) set(0,'DefaultAxesFontSize',18,'DefaultTextFontSize',18); subplot(2,2,1) plot(t,y(:,1),'r','linewidth',2); legend('B') subplot(2,2,2) plot(t,y(:,2),'r','linewidth',2); legend('G') subplot(2,2,3) plot(t,y(:,3),'r','linewidth',2); legend('C') subplot(2,2,4) plot(t,y(:,4),'r','linewidth',2); legend('N') function dydt=system(t,y) f=y(4)/(K_f+y(4)); g=y(3)/(K_g+y(3)); h=y(1)/(K_h+y(1));

    dydt1=-D*y(1)+mu_B*y(1)*min(f,g)-mu_G*h*y(2)-epsilon*y(1);
      dydt2=-D*y(2)+alpha_mo*mu_G*h*y(2)-w*y(2);
      dydt3=D*(C_in-y(3))-1/r*mu_B*min(f,g)+w*y(2)+epsilon*y(1);
       dydt4=D*(N_in-y(4))+theta_G*w*y(2)+(theta_B-theta_G*alpha_mo)*mu_G*h*y(2)-theta_B*epsilon*y(1)*min(f,g)+theta_B*mu_B*y(1);
        dydt=[dydt1;dydt2;dydt3;dydt4];
   end

end