How can I store values in empty matrix to plot later?

I am trying to simulate an epidemic model based on the SIR model. I am having trouble with my code. My graphs are not showing up like i think they should. It is not an equation error. I think there may be a problem with the storage of my values or something. Please help.

%Epidemic Simulation
%Author: James Metz
%Date: Nov 15, 2020

%Define Natural Paramteres:
p = 0.0144; %Natural Birth Rate (Davidson County)
u = 0.008; %Natural Death Rate (Davidson County)
f = 0.03; %Infection Rate
r = 0.075; %Recovery Rate
m = 0.0001; %Death due to Infection Rate
v = 0.092; %Vaccination Rate

%Define Evaluation Time Paramteres:
dt = 1; %Time increments (days)
tEnd = 365; %Simulation length (days)
t = 1:dt:tEnd;

%Initialize Variables: 
I = 10; %Initial infected population
R = 0; %Initial recovered population
S = 692587; %Initial susceptible population
I_1d = zeros(1, tEnd); 
R_1d = zeros(1, tEnd); 
S_1d = zeros(1, tEnd); 

S_1d(1) = 692587;
I_1d(1) = 10;
R_1d(1) = 0; 

%Loop through times:
for idx = 1:dt:tEnd-1
    %Initialization:
    S = S_1d(idx); 
    R = R_1d(idx); 
    I = I_1d(idx); 
    %Find changes in population numbers
    dS_dt = -f*S*I - S*u - S*v + S*p; %Drop in unifected population
    dI_dt = f*S*I - r*I - m*I - u*I; %Drop in infected population
    dR_dt = r*I + S*v - u*R; %Gain in Recovered population
    %Store new values:
    S_1d(idx+1) = S_1d(idx) + dS_dt; 
    R_1d(idx+1) = R_1d(idx) + dR_dt; 
    I_1d(idx+1) = I_1d(idx) + dI_dt; 
end

figure(1)
subplot(2,2,1)
plot(t, S_1d)
xlabel('Time (days)')
ylabel('Susceptible Population')
title('Drop in Susceptible Population due to Infection')
subplot(2,2,2)
plot(t, R_1d)
xlabel('Time (days)')
ylabel('Recovered Population')
title('Recovery Rate of Infected Persons')
subplot(2,2,3)
plot(t, I_1d)
xlabel('Time (days)')
ylabel('Infected Population')
title('Infection Rate of Susceptible Persons')
subplot(2, 2, 4)
plot(t, S_1d, 'green')
plot(t, R_1d, 'Blue')
plot(t, I_1d, 'red')

This is how my graphs are turining out:

[KSSV]

Check the values......are you getting NaN's after some point of the loop.

[James Metz]

Used isnan() function, but nothing yeilded a 'true'

Answer 1

[Alan Stevens]

You need to wotk with normalized population numbers. You can renormalize at the end.

%Epidemic Simulation
%Author: James Metz
%Date: Nov 15, 2020

%Define Natural Paramteres:
p = 0.0144; %Natural Birth Rate (Davidson County)
u = 0.008; %Natural Death Rate (Davidson County)
f = 0.03; %Infection Rate
r = 0.075; %Recovery Rate
m = 0.0001; %Death due to Infection Rate
v = 0.092; %Vaccination Rate

%Define Evaluation Time Paramteres:
dt = 1; %Time increments (days)
%tEnd = 365; %Simulation length (days)
tEnd =365;
t = 0:dt:tEnd;
elems = numel(t);

%Initialize Variables: 
I0 = 10; %Initial infected population
R0 = 0; %Initial recovered population
S0 = 692587; %Initial susceptible population
I_1d = zeros(1, elems); 
R_1d = zeros(1, elems); 
S_1d = zeros(1, elems); 

S_1d(1) = 1;              % Normalize the numbers
I_1d(1) = I0/S0;
R_1d(1) = R0/S0; 

%Loop through times:
for idx = 1:elems-1
    %Initialization:
    S = S_1d(idx); 
    R = R_1d(idx); 
    I = I_1d(idx); 
    %Find changes in population numbers
    dS_dt = -f*S*I - S*u - S*v + S*p; %Drop in unifected population
    dI_dt = f*S*I - r*I - m*I - u*I; %Drop in infected population
    dR_dt = r*I + S*v - u*R; %Gain in Recovered population
    %Store new values:
    S_1d(idx+1) = S + dS_dt*dt; 
    R_1d(idx+1) = R + dR_dt*dt; 
    I_1d(idx+1) = I + dI_dt*dt; 
end
S_1d = S_1d*S0;     % Renormalize the numbers
R_1d = R_1d*S0;
I_1d = I_1d*S0;
figure(1)
subplot(2,2,1)
plot(t, S_1d)
xlabel('Time (days)')
ylabel('Susceptible Population')
title('Drop in Susceptible Population due to Infection')
subplot(2,2,2)
plot(t, R_1d)
xlabel('Time (days)')
ylabel('Recovered Population')
title('Recovery Rate of Infected Persons')
subplot(2,2,3)
plot(t, I_1d)
xlabel('Time (days)')
ylabel('Infected Population')
title('Infection Rate of Susceptible Persons')
subplot(2, 2, 4)
plot(t, S_1d, 'green')
plot(t, R_1d, 'Blue')
plot(t, I_1d, 'red')

This results in

[James Metz]

I'm sorry, I don't understand what you did by normalizing the number? What does that mean/how do you do it?

[Alan Stevens]

Normalizing the population numbers means dividing the values of S, I and R by the total number of people in the population. Then your equations work on fractions of the population. So a fraction of 1 means the whole population, a fraction of 0.9 means 90% of the population and so on.

[James Metz]

Oh that makes sense. If you don't mind me asking, why is it necessary in a code?

[Alan Stevens]

It wouldn't be necessary if your equations were linear, but terms like S*I make it so.