In[]:=
ClearAll["`*"];Remove["`*"];
In[]:=
δ=0;​​γ=
1
10
;​​ω=1;​​s=ParametricNDSolve[{Subscript[ρ,22]'[t]==-γ*Subscript[ρ,22][t]-(I/2*ω*Subscript[ρ,21][t]-I/2*Conjugate[ω]*Subscript[ρ,12][t]),Subscript[ρ,11]'[t]==γ*Subscript[ρ,22][t]+(I/2*ω*Subscript[ρ,21][t]-I/2*Conjugate[ω]*Subscript[ρ,12][t]),Subscript[ρ,21]'[t]==-(I*δ+g)*Subscript[ρ,21][t]-I/2*ω*(Subscript[ρ,22][t]-Subscript[ρ,11][t]),Subscript[ρ,21]'[t]==Conjugate[Subscript[ρ,12]'[t]],Subscript[ρ,11][0]==1,Subscript[ρ,22][0]==0,Subscript[ρ,12][0]==0,Subscript[ρ,21][0]==0},{Subscript[ρ,11],Subscript[ρ,12],Subscript[ρ,21],Subscript[ρ,22]},{t,0,100},{g}];​​Plot[Evaluate[Table[Subscript[ρ,22][g][t]/.s,{g,0,2,0.3}]],{t,0,100},PlotRange->All]
Out[]=
20
40
60
80
100
0.2
0.4
0.6
0.8