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