giant spikes in model from equilibrium when run for relatively long time periods

Hi, I'm working on a model of Daphnia and algae in a chemostat, and when I run the code below for more than 300 days or so, the model will reach equilibrium within 100-200 days but then there are frequently giant spikes around 200-300 days that shoot up very quickly and then return to the equilibrium value. I showed this to my advisor and he thought they may be the result of numerical errors. I've tried decreasing the relative tolerance and max step (especially since my parameters range from 1 to 1.28e-8), but if I decrease the max step more than I have here (1e-2) it takes too long to run (i.e. hasn't produced a graph even if left overnight). Any suggestions on ways to fix this would be much appreciated! Thank you!

function algae_daphnia_NOrecycling()

% Time Variable T = [0 500];

% Solver Options options = odeset('RelTol', 1E-12, 'MaxStep', 1E-2);

%state variables N_Ainit = 0.1; Q_Ainit = 0.00000027778; N_dinit = 0.1420; R_init = 2.5E-6; x = [N_Ainit; Q_Ainit; N_dinit; R_init];

% Solver [T x] = ode23s(@odefun_dynamic, T, x, options);

%Plot results - four graphs figure() subplot(2,2,1), plot(T, x(:,1)); xlabel('time (days)'); ylabel('Algae biomass (mg-C/L)'); subplot(2,2,2), plot(T, x(:,2)); xlabel('time (days)'); ylabel('Algae quota (Mol-P/mg-C)'); subplot(2,2,3), plot(T, x(:,3)); xlabel('time (days)'); ylabel('Daphnia biomass (mg-C/L)'); subplot(2,2,4), plot(T, x(:,4)); xlabel('time (days)'); ylabel('Available resource (Mol-P/L)');

end

function x_prime = odefun_dynamic(T, x)

% flow rate parameter F = 0.1;

% set parameters S = 2.5e-6; %initial substrate concentration p_max = 3.38e-6; %max intake rate of nutrient (P) by algae K_p = 1.29e-8; %half saturation constant of nutrient intake by algae u_Amax = 1; %apparent max growth rate Q_Amin = (2.5E-7); %subsistence quote for algae q_max = 1; %max uptake rate of algae by Daphnia K_q = 0.164; %half saturation constsant of uptake of algae by Daphnia gamma = 0.5; %fraction of biomass of algae assimilated by Daphnia m_d = 0.02; %mortality of Daphnia

% Get state variables from x N_A = x(1); Q_A = x(2); N_d = x(3); R = x(4);

%functions

M_A = F * N_A;

p = (p_max * R) / (K_p + R);

I_R = N_A * p;

r_A = (u_Amax * N_A) * ((Q_A - Q_Amin) / (Q_A));

M_d = m_d * N_d;

I_A = (N_d * q_max * N_A) / (K_q + N_A);

r_d = gamma * I_A;

% differential equations

N_A_prime = r_A - M_A - I_A;

Q_A_prime = p - ((Q_A - Q_Amin) * u_Amax);

N_d_prime = r_d - M_d;

R_prime = (F * S) - (F * R) - I_R;

% Return derivatives in a single column vector x_prime = [N_A_prime; Q_A_prime; N_d_prime; R_prime];

end