clear all;
close all;
% Parameters of the system:
N = 10000;
a = 6.90675*10^(-5);
b = 10;
r = 0.1;
% Initial data of the system:
t0 = 0;
I0 = 20;
D0 = 0;
M0 = 0;
S0 = N - I0;
% Time interval and step size for a numerical solution:
T = 50;
h = 0.1;
n = T/h;
% Iterative loop to define a numerical solution in discrete times:
t(1) = 0; S(1) = S0; D(1) = D0; M(1) = M0; I(1) = I0;
for k = 1 : n
t(k+1) = t(k) + h;
S(k+1) = S(k) - h*a*S(k)*I(k);
D(k+1) = D(k) + h*r*I(k);
M(k+1) = M(k) + h*b;
I(k+1) = I(k) + h*a*S(k)*I(k) - h*r*I(k) - h*b;
end
% Showing results as functions of time (1)-(4) and on the phase plane (5)
figure(1);
plot(t,S); xlabel('days'); ylabel('Number of Susceptible people');
figure(2);
plot(t,I); xlabel('days'); ylabel('Number of Infected people');
figure(3);
plot(t,M); xlabel('days'); ylabel('Number of medicated people');
figure(4);
plot(t,D); xlabel('days'); ylabel('Number of dead people');
figure(5);
plot(I,S); xlabel('Number of Infected people'); ylabel('Number of Susceptible people');
% Finding the number of survived people after dT = 15 days
dT = 15;
nn = dT/h;
R = N - D(nn+1);
fprintf('The number of survived people after %4.0f days is %4.0f\n',t(nn+1),R);