The simulation shows some results, although I didn't check/compare your equations with the standard models in
Note that ϵ is not given. If 
, then 
and 
converge very fast.
, then and converge very fast.function Demo_SIR
close all;
clear;
clc;
Tspan = [0 1e3]; % simulation time
y0 = [1000; 500; 5; 5; 0; 0]; % initial values
[t, y] = ode45(@SIR, Tspan, y0);
plot(t, y, 'linewidth', 1.5)
grid on
xlabel('Time')
ylabel('Individual')
legend({'$S_{1}$', '$S_{2}$', '$I_{1}$', '$I_{2}$', '$R_{1}$', '$R_{2}$'}, 'Interpreter', 'latex', 'location', 'best')
title('Responses of the SIR model')
end
function dydt = SIR(t, y)
dydt = zeros(6,1);
alpha = 0;
beta = 0.00002;
gamma = 0.01;
epsilon = 0.125;
V = 0.4;
S1 = y(1);
S2 = y(2);
I1 = y(3);
I2 = y(4);
R1 = y(5);
R2 = y(6);
dydt(1) = - beta*S1*((1 - epsilon)*(1 + alpha/V)*I1 + (1 - epsilon)*I2);
dydt(2) = - beta*S2*(I1 + I2);
dydt(3) = - (gamma*I1) - (- beta*S1*((1 - epsilon)*(1 + alpha/V)*I1 + (1 - epsilon)*I2));
dydt(4) = - (gamma*I2) - (- beta*S2*(I1 + I2));
dydt(5) = gamma*I1;
dydt(6) = gamma*I2;
end
Result:
