clc;clear;close?all%?計算意大利
hubei_data=xlsread（‘意大利數據.xlsx‘）;
hubei_data=hubei_data（13:end:）;?%過濾2月14號以前數據
%%?參數
Q=1;????%?潛伏者相對于感染者傳播能力的比值
A=1/14;??%隔離時間14天
G=1/7?;?%?潛伏期向感染者的轉化速率潛伏期為7d
c=3.2;??%接觸率
p=1;???%有效接觸系數
Si=0.13;??%感染者隔離速率
Sq=0.13;??%?隔離
YI=0.007;YH=0.014;?B=2.05*10^（-9）;
q=10^（-6）;a=0.0399;
h=0.01?;??%求解的步長
T=60?;??%預測50天的
t=0:h:T;
S=60430000?;E=3;I=3*1.5;SQ=9;EQ=2;H=I+EQ;R=42;
DS=zeros（1length（t））;?????DS1=zeros（1length（t））;?????%?參數保存地方
DE=zeros（1length（t））;???DE1=zeros（1length（t））;
DI=zeros（1length（t））;???DI1=zeros（1length（t））;
DSQ=zeros（1length（t））;??DSQ1=zeros（1length（t））;
DEQ=zeros（1length（t））;???DEQ1=zeros（1length（t））;
DH=zeros（1length（t））;???DH1=zeros（1length（t））;
DR=zeros（1length（t））;???DR1=zeros（1length（t））;
DS（1）=S;%參數初始化
DE（1）=E;DI（1）=I;DSQ（1）=SQ;DEQ（1）=EQ;DH（1）=H;DR（1）=R;
for?i=2:length（t）
????DS1（i-1）=-（p*c*B+p*c*q*（1-B））*（?DS（i-1）?）*（?DI（i-1）+1*DE（i-1）??）+A*DSQ（i-1）;
????DS（i）=DS（i-1）+DS1（i-1）*h;
????DE1（i-1）=p*c*B*（1-q）*DS（i-1）*（?DI（i-1）+1*DE（i-1）?）-G*DE（i-1）;
????DE（i）=DE（i-1）+DE1（i-1）*h;
????DI1（i-1）=G*DE（i-1）-（Si+a+YI）*DI（i-1）;
????DI（i）=DI（i-1）+DI1（i-1）*h;
????DSQ1（i-1）=p*c*q*（1-B）*DS（i-1）*（?DI（i-1）+1*DE（i-1）?）-A*DSQ（i-1）;
????DSQ（i）=DSQ（i-1）+DSQ1（i-1）*h;
????DEQ1（i-1）=p*c*B*q*DS（i-1）*（?DI（i-1）+1*DE（i-1）?）-G*DE（i-1）;
????DEQ（i）=?DEQ（i-1）+?DEQ1（i-1）*h;
????DH1（i-1）=Si*DI（i-1）+Sq*DEQ（i-1）-（a+YH）*DH（i-1）;
????DH（i）=DH（i-1）+DH1（i-1）*h;
????DR1（i-1）=（YI）*DI（i-1）+YH*DH（i-1）;
????DR（i）=DR（i-1）+DR1（i-1）*h;
end

figure（1）;plot（tDI）;?hold?on;plot（t-DR）;?hold?on?;?
plot（hubei_data（:1）-12hubei_data（:3）‘*‘）;?hold?on;
plot（hubei_data（:1）-12hubei_data（:4）‘*‘）;?hold?on;
xlabel（‘時間-初始時間是2月14日‘）;ylabel（‘人數（人）‘）;
legend（‘當前實際感染人數‘‘當前實際治愈人數‘‘理論預測-當前感染‘‘理論預測-當前治愈‘）;
title（‘意大利新冠肺炎感染趨勢‘）;grid?on