% Parameters
c = 0.0020;
a = 0.02;
% SIR model (Ordinary Differential Equations)
odefcn = @(t,x) [- c*x(1)*x(2);
c*x(1)*x(2) - a*x(2);
a*x(2)];
%%The t in the defintion of odefcn feels redundant
% Simulation time span
t0 = 0;
h = 0.1; % step size
tf = 100;
t = t0:h:tf;
% Initial values
y0 = [100 1 1];
% Calling RK4 Solver is similar to ode45 solver
yy = rk4(odefcn, t, y0);
%Define array for legend
str = {'S(t)','I(t)','R(t)'};
%loop to plot in separate figures
for k=1:size(yy,1)
figure(k)
% Plotting the solutions
plot(t, yy(k,:), 'linewidth', 1.5)
%Change other properties accordingly
grid on
xlabel('Time in Days')
ylabel('Population')
title('Time responses of SIR model')
legend(str{k}, 'location', 'best')
end
% Code for Runge-Kutta 4th-order Solver
function y = rk4(f, x, y0)
y(:, 1) = y0; % initial condition
h = x(2) - x(1); % step size
n = length(x); % number of steps
for j = 1 : n-1
k1 = f(x(j), y(:, j)) ;
k2 = f(x(j) + h/2, y(:, j) + h/2*k1) ;
k3 = f(x(j) + h/2, y(:, j) + h/2*k2) ;
k4 = f(x(j) + h, y(:, j) + h*k3) ;
y(:, j+1) = y(:, j) + h/6*(k1 + 2*k2 + 2*k3 + k4) ;
end
end
