What is wrong with my SIR Model?

I tried to use SIR Model for my homework but it doesn't run. This is the question.

• Let S(t) denote the number of susceptible at time t (in hours).

• Let I(t) denote the number of infected individuals at time t (in hours).

• Let R(t) denote the number of recovered individuals at time t (in hours).

Use MATLAB’s or R’s built-in functions (i.e. ode45 or deSolve) to solve the following system of ODEs from t = 0 to t = 720 hours:

A) The spread of the disease is modeled as

dS/dt= -βSI,

dI/dt = βSI-γI,

dR/dt = γI,

with the initial conditions S(0) = 50, I(0) = 1, R(0) = 0 and γ = 0.0083, β = 0.0006.

clc;
clear;

function dPop=Diff_2_1(t, pop, parameter);

    beta=0.0006;
    gamma=0.0083;
    S=pop(1);
    I=pop(2);
    R=pop(3);
    
    dPop=zeros(3,1);

    dPop(1)= -beta*S*I;
    dPop(2)= beta*S*I - gamma*I;
    dPop(3)= gamma*I;
end
function [t,S,I,R] = Program_2_1(beta,gamma,S0,I0,MaxTime)
    S=S0;
    I=I0;
    R=1-S-I;
    [t, pop]=ode45(@(t,y) Diff_2_1(t,y,[beta gamma]),[0 MaxTime],[S I R]);
    S=pop(:,1); I=pop(:,2); R=pop(:,3);

[t,S,I,R] = Program_2_1(0.0006,0.0083,50,1,720);
plot(t,S,"-r",t,I,"-g",t,R,"-b")
xlim([0 200])

xlabel("Time","fontweight","bold")
ylabel("Number","fontweight","bold")
h = legend("S","I","R");
legend(h,"show")
end

The code doesn't run. I mean there is no output. Maybe the problem is because I ended the function in the wrong place but I still don't figure out it. Can you show me where is my mistake?

Answer 1

[praguna manvi]

Hi @Eymen Dogrusever,

As I see you have ended the second function wrongly at the function call "Program_2_1(0.0006,0.0083,50,1,720);" and also "legend" is used incorrectly, here is the corrected version of the above snippet :

clc;
clear;

function dPop = Diff_2_1(t, pop, parameter)
    beta = parameter(1);
    gamma = parameter(2);
    S = pop(1);
    I = pop(2);
    R = pop(3);
    
    dPop = zeros(3,1);
    dPop(1) = -beta * S * I;
    dPop(2) = beta * S * I - gamma * I;
    dPop(3) = gamma * I;
end

function [t, S, I, R] = Program_2_1(beta, gamma, S0, I0, MaxTime)
    R0 = 0; % Initial condition for R
    [t, pop] = ode45(@(t, y) Diff_2_1(t, y, [beta gamma]), [0 MaxTime], [S0 I0 R0]);
    S = pop(:,1);
    I = pop(:,2);
    R = pop(:,3);
    
    plot(t, S, '-r', t, I, '-g', t, R, '-b')
    xlim([0 200])
    xlabel('Time', 'fontweight', 'bold')
    ylabel('Number', 'fontweight', 'bold')
    legend('S', 'I', 'R') % Directly pass the strings for the legend
end

% Call the function
[t, S, I, R] = Program_2_1(0.0006, 0.0083, 50, 1, 720);

●