clear all;
close all;

% Initial data of the system:
N = 10000;
t0 = 0;
I0 = 20;
D0 = 0;
M0 = 0;
S0 = N - I0;
y0 = [S0,D0,M0,I0];

% Time interval and step size for a numerical solution:
T = 50;
h = 0.1; 
n = T/h; 

% Call to MATLAB ODE solver "ode23"
[t,y] = ode23(@epidem,[t0,T],y0);

% Showing results as functions of time (1)-(4) and on the phase plane (5)
figure(1); 
plot(t,y(:,1)); xlabel('days'); ylabel('Number of Susceptible people'); 
figure(2); 
plot(t,y(:,4)); xlabel('days'); ylabel('Number of Infected people'); 
figure(3); 
plot(t,y(:,3)); xlabel('days'); ylabel('Number of medicated people'); 
figure(4); 
plot(t,y(:,2)); xlabel('days'); ylabel('Number of dead people'); 
figure(5); 
plot(y(:,4),y(:,1)); xlabel('Number of Infected people'); ylabel('Number of Susceptible people'); 

% Finding the number of survived people after dT = 15 days
dT = 15;
[t,y] = ode23(@epidem,[t0,dT],y0);
R = N - y(length(t),2); 
fprintf('The number of survived people after %4.0f days is %4.0f\n',t(length(t)),R);
