MATLAB: I am trying to solve and simulate an epidemic model using ODE 45 with a guide provided from a textbook: Below are the list of functions i have defined:

modelode45

A = 13.69863014; B = 8.95890411; C = 1.585205479; D = 0.0095; beta = 0.000005937 F = 0.000095; G = 0.21; H =0.1; J = 0.0769230; K = 0.21; gamma = 0.11572; varpi = 0.771466; psi = 0.2; tau =0.00004466;
function dydt = SMVI(t,y)
dydt(1) = ( (1 – B – C) * A + varpi * y(2) + K * y(3) – (psi + gamma + D) * y(1) – beta * y(1) * y(4));
dydt(2) = ( (A * B ) + psi * y(1) – (varpi + J + D) * y(2) – G * beta * y(2) * y(4)) ;
dydt(3) = ((A * C) + gamma * y(1) + J * y(2) – (omega + D) * y(3) – H * beta * y(3) * y(4));
dydt(4) = (beta * y(1) * y(4) + G * beta * y(2) * y(4)… + H * beta * y(3) * y(4)… – (tau + F + D) *y(4));
dydt = [dydt(1);dydt(2);dydt(3);dydt(4)];
My problem now is how to create an m file to plot my graphs for each compartment as I do not seems to understand the text anymore.
All suggestions will be greatly appreciated.

Best Answer

You need to save your ‘SMVI’ function as SMVI.m on your MATLAB search path:
function dydt = SMVI(t,y)
    A = 13.69863014; 
    B = 8.95890411;    
    C = 1.585205479;      
    D = 0.0095;           
    beta = 0.000005937;    
    F = 0.000095;      
    G = 0.21;        
    H =0.1;      
    J = 0.0769230;    
    K = 0.21;         
    gamma = 0.11572;  
    varpi = 0.771466; 
    psi = 0.2;
    tau =0.00004466;
    omega = 0.42;                           % <— PROVIDE CORRECT VALUE
dydt(1) = (  (1 - B - C) * A + varpi * y(2) + K * y(3) - (psi + gamma + D) * y(1) - beta * y(1) * y(4));            
dydt(2) = ( (A * B ) + psi * y(1) - (varpi + J + D) * y(2) - G * beta * y(2) * y(4));            
dydt(3) = ((A * C) + gamma * y(1) + J * y(2) - (omega + D) * y(3) - H * beta * y(3) * y(4));
dydt(4) = (beta * y(1) * y(4) + G * beta * y(2) * y(4)... 
        + H * beta * y(3) * y(4)... 
        - (tau + F + D) *y(4));
dydt = [dydt(1);dydt(2);dydt(3);dydt(4)];
end
Then create a script file to integrate it, for example as:
tspan = [0 100];
y0 = ones(1,4) * 1E-8;                      % <— PROVIDE CORRECT INITIAL CONDITIONS
[T,Y] = ode45(@SMVI, tspan, y0);
figure(1)
plot(T, Y)
grid
legend('C_1', 'C_2', 'C_3', 'C_4', 'Location','SE')
Experiment to get the results you want.
Related Question