I would appreciate if someone can help with the following issue. I have the following ODE:
this is the project that I have worked to create a Matlab code to it
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7376536/#Equ1
the human population N(t)=S(t)+E(t)+IA(t)+IS(t)+R(t).
Model parameter name Symbol Value
Birth rate of the human population b 0.00018 days−1
Natural human death rate μ 4.563×10−5 days−1
Human life expectancy 1μ 21915 days or 60 years
Natural death rate of pathogens in the environment μP 0.1724 days−1
Life expectancy of pathogens in the environment 1μP 5.8 days
Proportion of interaction with an infectious environment α1 0.10
Proportion of interaction with an infectious individual α1 0.10
Rate of transmission from S to E due to contact with P β1 0.00414
Rate of transmission from S to E due to contact with IA and/or IS β2 0.0115
Proportion of symptomatic infectious people δ 0.7
Progression rate from E back to S due to robust immune system ψ 0.0051
Progression rate from E to either IA or IS ω 0.09
Death rate due to the coronavirus σ 0.0018
Rate of recovery of the symptomatic population γS 0.05 days−1 or 120days
Rate of recovery of the asymptomatic human population γA 0.0714 days−1
Rate of virus spread to environment by symptomatic infectious individuals ηS 0.1 days−1 or 110days
Rate of virus spread to environment by asymptomatic infectious individuals ηA 0.05 days−1 or 120days
I'm trying to solve a system of coupled ODEs using a 4th-order Runge-Kutta method for my project work in Matlab .
I'm stuck in first loop , can't understand what next steps to use k1,k2,k3,k4 in the method
function y = code(S0,E0,IA0,IS0,R0,P0)
% Birth rate of the human population
b= 0.00018;
% Natural human death rate
u = 0.00004563
% Human Life expectancy
inverse_u = 1/u;
% Natural deapth rate of pathogens in the environment
up= 0.1724;
%Life expectancy of pathogens in environment
inverse_up = 1/up;
%Proportion of interaction with an infectious environment
a1= 0.1;
% Proportion of interaction with an infectious individual
a2= 0.1;
% Rate of transmission from S to E due to contact with P
B1= 0.00414;
% Rate of transmission from S to E due to contact with Ι_Α and/or Ι_S
B2= 0.0115;
% δ Proportion of symptomatic infectious people
Delta = 0.7;
%ψ Progression rate from E back to S due to robust immune system
psi = 0.0051;
% ω Progression rate from E to either Ι_Αor Ι_S
Omega = 0.09;
% σ Death rate due to the Coronavirus
Sigma = 0.0018
% νs Rate of recovery of the symptomatic population
Vs = 0.05;
% νa Rate of recovery of the symptomatic human population
Va = 0.0714;
% ηs Rate of virus spread to environment by symptomatic infectious individuals
Eta_s = 0.1;
% ηa Rate of virus spread to environment by asymptomatic infectious individuals
Eta_a = 0.1;
test = -1;
delta = 0.001;
M = 90;
t=linspace(0,10,M+1);
h=1/M;
h2 = h/2;
S=zeros(1,M+1);
E=zeros(1,M+1);
IA=zeros(1,M+1);
IS=zeros(1,M+1);
R=zeros(1,M+1);
P=zeros(1,M+1);
S(1)=S0;
E(1)=E0;
IA(1)=IA0;
IS(1)= IS0;
R(1)=R0;
P(1)=P0;
%() = () + () +() + () + ().
N(1)=S0+E0+IA0+IS0+R0
lambda1=zeros(1,M+1);
lambda2=zeros(1,M+1);
lambda3=zeros(1,M+1);
lambda4=zeros(1,M+1);
lambda5=zeros(1,M+1);
while(test < 0)
oldP = P;
oldS = S;
oldE = E;
oldIA = IA;
oldIS = IS;
oldR = R;
oldN = N;
oldlambda1 = lambda1;
oldlambda2 = lambda2;
oldlambda3 = lambda3;
oldlambda4 = lambda4;
oldlambda5 =lambda5;
for i = 1:M
m11 = b -( ( B1 *S(i)*P(i) ) / (1+a1*P(i) ) ) - ( ( B2*S(i) *(IA(i) + IS(i) ) ) / ( 1+a2*(IA(i) + IS(i)) ) ) + psi*E(i) - u*S(i) ;
m12 = ( ( B1 *S(i)*P(i) ) / (1+a1*P(i) ) ) + ( ( B2*S(i) *(IA(i) + IS(i) ) ) / ( 1+a2*(IA(i) + IS(i)) ) ) -psi*E(i) - u*E(i) - Omega *E(i);
m13 = (1 - Sigma ) *Omega*E(i) - (u + Sigma ) * IA(i) - Va *IA(i) ;
m14 = Sigma * Omega*E(i) - ( u + Sigma ) * IS(i) - Vs*IS(i) ;
m15 = Vs*IS(i) + Va *IA(i) - u*R(i) ;
m16 = Eta_a * IA(i) + Eta_s * IS(i) - up*P(i);
