My answer

In[]:=
Clear["Global`*"]
In[]:=
(*Recoveryrate*)​​γ=0.1;​​(*infectionrate*)​​β=0.4;​​tMax=90;​​n=10^6;​​​​(*FoI*)​​λ[t_]=βi[t]/n;​​(*Definethesystemofdifferentialequations*)​​eqns={s'[t]==-λ[t]*s[t],i'[t]==λ[t]*s[t]-γ*i[t],r'[t]==γ*i[t],s[0]==999999,i[0]==1,r[0]==0};​​(*Solvethesystem*)​​sol=NDSolve[eqns,{s,i,r},{t,0,tMax}];​​(*Plotthesolutions*)​​Plot[Evaluate[{s[t],i[t],r[t]}/.sol],{t,0,tMax},PlotLegends->{"Susceptible","Infected","Recovered"},AxesLabel->{"Time","Population"},PlotStyle->{Blue,Red,Green}]
Out[]=
Susceptible
Infected
Recovered
In[]:=
Reff[t_]=(β/γ)*(s[t]/n)
Out[]=
4.×
-6
10
s[t]
In[]:=
Plot[Reff[t]/.sol,{t,0,90},PlotStyle->Gray]
Out[]=
In[]:=
inc[tmax_]:=NIntegrate[{λ[t]*s[t]}/.sol,{t,0,tmax}]
In[]:=
Plot[inc[t],{t,7,14},AxesLabel->{"Time","Incidence"},PlotStyle->Orange]
Out[]=
In[]:=
inc[14]-inc[7]/.sol
Out[]=
{{{78.018}}}

colleague answer

In[]:=
Clear["Global`*"];
In[]:=
​​γ=0.1;​​β=0.41;​​tMax=90;​​n=10^6;​​λ[t]:=βi[t]/n;​​eqns={s'[t]==-λ[t]s[t],i'[t]==λ[t]s[t]-γ*i[t],r'[t]==γ*i[t],incident'[t]==λ[t]*s[t],s[0]==999999,i[0]==1,r[0]==0,incident[0]==0};​​sol=NDSolve[eqns,{s,i,r,incident},{t,0,tMax}];​​Plot[Evaluate[{s[t],i[t],r[t],incident[t]}/.sol],{t,0,tMax},PlotLegends->{"Susceptible","Infected","Recovered","Incidence"},AxesLabel->{"Time","Population"},PlotStyle->{Blue,Red,Green,Directive@{Black,Dashed}}]
Out[]=
Susceptible
Infected
Recovered
Incidence
In[]:=
incident[14]-incident[7]/.sol​​
Out[]=
{89.8569}