seir分段函数代码

clear;clc;
N=17560000;
I=1;%传染者
R=0;%康复者
D=0;%死亡患者数量
E=0;%潜伏者
S=N-I;%易感染者
m=4;%接触病患的人数（一阶段为未防控，剩余阶段参数递减的有：病患接触人数m,
%                     潜伏者接触人数r和患病概率B.)
m1=1;%二阶段病患接触
m2=0;%三四阶段病患接触
a=0.125;%潜伏者患病概率
B=0.6;%感染概率
B00=0.3;
B01=0.2;
B02=0.1
y=0.143;%康复概率
B2=0.03;%接触潜伏者
B3=0.859;%转阴率
r2=8;%潜伏者每天接触的人数
r4=3;%2阶段潜伏者
r5=2;%三阶段
r6=1;%四阶段
k=0.025373;%死亡率
T=0:200;
for idx =1:length(T)-1
    if idx<30
S(idx+1)=S(idx)-m*B*I(idx)-r2*B2*E(idx)+B3*E(idx);%易感人数迭代  
E(idx+1)=E(idx)+m*B*I(idx)+r2*B2*E(idx)-B3*E(idx);%潜伏者人数迭代 
I(idx+1)=I(idx)+a*E(idx)-(y+k)*I(idx);%患病人数迭代              
R(idx+1)=R(idx)+y*I(idx);%康复人数迭代
D(idx+1)=D(idx)+k*I(idx);%死亡患者人数迭代
    elseif idx>=30&&idx<40
S(idx+1)=S(idx)-m1*B00*I(idx)-r4*B2*E(idx)+B3*E(idx);%易感人数迭代  
E(idx+1)=E(idx)+m1*B00*I(idx)+r4*B2*E(idx)-B3*E(idx);%潜伏者人数迭代 
I(idx+1)=I(idx)+a*E(idx)-(y+k)*I(idx);%患病人数迭代              
R(idx+1)=R(idx)+y*I(idx);%康复人数迭代
D(idx+1)=D(idx)+k*I(idx);%死亡患者人数迭代
    elseif idx>=40&&idx<=50
S(idx+1)=S(idx)-m1*B01*I(idx)-r5*B2*E(idx)+B3*E(idx);%易感人数迭代  
E(idx+1)=E(idx)+m1*B01*I(idx)+r5*B2*E(idx)-B3*E(idx);%潜伏者人数迭代 
I(idx+1)=I(idx)+a*E(idx)-(y+k)*I(idx);%患病人数迭代              
R(idx+1)=R(idx)+y*I(idx);%康复人数迭代
D(idx+1)=D(idx)+k*I(idx);%死亡患者人数迭代
     elseif idx>=50&&idx<=length(T)-1
S(idx+1)=S(idx)-m2*B02*I(idx)-r6*B2*E(idx)+B3*E(idx);%易感人数迭代  
E(idx+1)=E(idx)+m2*B02*I(idx)+r6*B2*E(idx)-B3*E(idx);%潜伏者人数迭代 
I(idx+1)=I(idx)+a*E(idx)-(y+k)*I(idx);%患病人数迭代              
R(idx+1)=R(idx)+y*I(idx);%康复人数迭代
D(idx+1)=D(idx)+k*I(idx);%死亡患者人数迭代 
    end
end
plot(T,E,'--r',T,I,'-k');
grid on;
xlabel('日期');
ylabel('人数/人');
title('SEIR  模型');
%plot(T,E,T,I,T,R,T,D);
plot(T,E,'--r',T,I,'-k');
grid on;
xlabel('时间/天');
ylabel('人数/人');
axis([0 100 0 200]);
%legend('潜伏者','感染者','康复者','死亡者');
legend('潜伏者','感染者');
set(gca,'ytick',[0:50:500]);
%title('疫情情况');