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);
