Вмешательство в оду решатель, использующий для цикла

Question

Я пытаюсь запустить систему уравнений для моделирования популяций пчел и клещей с инфекциями. Тем не менее, как только популяция пчел достигнет определенного порога размера, я хочу, чтобы половина пчел ушла, а также определенный процент клещей тоже. Я пытаюсь запустить цикл проверки размера популяции каждый раз. Мой код функции:

function dxdt = project_v3(t,x,par)

dxdt = zeros(4,1); %pre-allocation of the dependent variables

m = x(1); %infected mites
z = x(2); %healthy bees
y = x(3); %infected bees
M = x(4); %total mites

g1 = par(1); %gamma 1 in model
g2 = par(2); %gamma 2
n = par(3);

Mu = [500 1500 500 0]; %vectors for each function of t, with values per season starting with spring
B1 = [0.1593 0.1460 0.1489 0.04226]; %[spring summer fall winter]
B2 = [0.04959 0.03721 0.04750 0.008460];
B3 = [0.1984 0.1460 0.1900 0.03384];
d1 = [0.02272 0.04 0.02272 0.005263];
d2 = [0.2 0.2 0.2 0];
k = [0.000075 0.00003125 0.000075 0.000075];
r = [0.0165 0.0165 0.0045 0.0045];
alpha = [0.4784 0.5 0.5 0.4784]; %alpha
K = [8000 12000 8000 6000];

S = floor(mod(t,365)/91.25)+1; 

%System of ODEs
dxdt(1) = B1(S).*(M-m)*(y/(z+y))-B2(S).*m*(z/(z+y));
dxdt(2) = Mu(S).*(z^n/(K(S).^n+z^n))*(exp(1)^(-m.*k(S)))-B3(S).*m*(z/(z+y))-d1(S).*z-g1*M*z; %Use Kh for K high and Kl for K low
dxdt(3) = B3(S)*m*(z/(z+y))-d2(S)*y-g2*M*y;
dxdt(4) = r(S).*M*(1-(M/(alpha(S).*(z+y))));

Это хорошо работает, когда я решаю это. Теперь я пытаюсь использовать решатель с циклом for для корректировки размеров популяции, если они достигают определенного размера:

close all;
clear all;

Y = 3; %number of years to run for
Tmax = Y*365; %years times days
tspan = 1:1:Tmax; %tspan from day 1 to Tmax

Mu = [500 1500 500 0]; %vectors for each function of t, with values per season starting with spring
B1 = [0.1593 0.1460 0.1489 0.04226]; %[spring summer fall winter]
B2 = [0.04959 0.03721 0.04750 0.008460];
B3 = [0.1984 0.1460 0.1900 0.03384];
d1 = [0.02272 0.04 0.02272 0.005263];
d2 = [0.2 0.2 0.2 0.005300];
k = [0.000075 0.00003125 0.000075 0];
r = [0.0165 0.0165 0.0045 0.0045];
alpha = [0.4784 0.5 0.5 0.4784];
K = [8000 12000 8000 6000];

Day = zeros(Tmax);
count = 1;

for t = 1:1:Tmax
S = floor(mod(t,365)/91.25)+1; %Denotes the season
Day(t) = Mu(S);
end

g1 = .0000001 ; %gamma 1 in model
g2 = .0000001; %gamma 2 in model
n = 2; % n >1; we assume n = 2

par = [g1 g2 n];

% Define the Initial conditions
m0 = 0;
z0 = 20000;
y0 = 0;
M0 = 100;

x0 = [m0 z0 y0 M0]; % Creates vector of initial conditions

for i = 1:1:Tmax
    c = .95; %threshold of maximum
    Zmax = 20342; %Max Z at disease-free equilibrium    
    a = .35; %proportion of mites staying
    b = .5; %proportion of uninfected bees staying  

    options = odeset('RelTol', 1e-4, 'NonNegative', [1 2 3 4]); %Converts all negative populations to 0
    [t,x] = ode45(@project_v3,tspan,x0,options,par);

    if x(i,2) >= c*Zmax
       x(i,2)=b*x(i,2);
       x(i,4)=a*x(i,4);
       x(i,1)=a*x(i,1);
    end

figure(1)

h1=plot(t,x);
xlabel('Time in Days')
ylabel('Number of Individuals')
title('Honeybees & Varroa Mites')
legend('Infected Mites','Healthy Bees','Infected Bees','Total Mites')
set(gca,'FontSize',10,'Linewidth',1.5)
set(h1,'Linewidth',1.5)
end

Тем не менее, он не строит графики с настроенными популяциями в определенное время, а скорее динамически пробегает предварительно построенный график и корректирует размеры по мере их поступления. Любой совет? Спасибо!