ode45 solver SIR model

I've used ode45 to solve a simple SIR model, I've got the graph to work as I wish but I'm strugling to output any numerical values to discuss.

I'm wanting to find where dI/dt = 0, for the time where the pandemic will be at its peak and the area under the curve for the total number of infected.

Here's my (probably very messy!) code:

%user parameters
N = 45742000; %total population
I0 = 1; %initial infected population
r0 = 12.9 %reproductive value
R0 = 0;%initial recovered population
i_period = 9; %duration of infectious period
tspan = [1, 50]; %length of time to run simulation over


%rate constant calculations
mu = 1 / i_period %recovery rate
beta = r0 * mu / N %infection rate

S0 = N - I0 - R0 %initial susceptible population

N0 = I0 + R0 + S0; %total population

%---------------------------------------------------
%feeding parameters into function
pars = [beta, mu, N, r0];
y0 = [S0 I0 R0];

Running SIR model function 
%using the ode45 function to perform intergration 
[t,y] = ode45(@sir_rhs, tspan, y0, [], pars);
figure()
plot(t,y(:,2), 'r');
xlim(tspan);
ylabel('Population (n)');
xlabel('Time (days)');
legend('Infected','Location','SouthEast');

function f = sir_rhs(t,y,pars)
f = zeros(3,1);
f(1) = -pars(1)*y(1)*y(2);
f(2) = pars(1)*y(1)*y(2) - pars(2)*y(2);
f(3) = pars(2) * y(2);
end

---------------------------------------------------------------------------------------------------

function f = sir_rhs(t, y, pars); %function contains differential equations
f = zeroes(4, 1); %creates an empty matrix which will be filled with values for susceptible, infected, recovered and total population, respectively 
beta = pars(1); 
mu = pars(2); 
N = pars(3); 
r = pars(4); 
S= y(1); 
I= y(2); 
R= y(3); 
n = y(4); 
f(1) = -beta*S*I; %susceptible population differential equation
f(2) = beta*S*I-mu*I; %infected populatin differential equation
f(3) = mu * I;  %recovered population differential equation
end

Im only interested in solving f(2) where it's equal to zero, which I think will be where the disease is at maximum or where it ends (I don't think it reaches zero more than once here because there's a single peak. I could use the sum of y(2) at each time point to calculate the total infected throughout the time course but there must be a simpler alternative!

Answer 1

[Mischa Kim]

Hi Ollie,

to get the area under the curve, just add another integration variable that integrates over the curve, y(2), I called it f(4). To locate the peak of the curve use an events function as shown below. ve contains the time stamp of the peak.

y0 = [S0 I0 R0 0];
%using the ode45 function to perform intergration 
options = odeset('Events',@events);
[t,y,ve,~,~] = ode45(@sir_rhs, tspan, y0, options, pars);
figure()
plot(t,y(:,2),t,y(:,4));
xlim(tspan);
ylabel('Population (n)');
xlabel('Time (days)');
legend('Infected','Location','SouthEast');

function f = sir_rhs(~,y,pars)
f    = zeros(4,1);
f(1) = -pars(1)*y(1)*y(2);
f(2) = pars(1)*y(1)*y(2) - pars(2)*y(2);
f(3) = pars(2) * y(2);
f(4) = y(2);
end

function [value,isterminal,direction] = events(~,y,pars)
value      = pars(1)*y(1)*y(2) - pars(2)*y(2); 
isterminal = 0;   % Do not stop the integration when zero is detected
direction  = 0;   % Detect all zeros
end