Numerical Solutions of SIT in continuous time; here, evolution of a preference; haploid model where nmi = nfi
and sterile releases are proprotional

Pars={(*S -> 11.745*^6*)S-> 3*^6,b-> 10,K-> 1.0*^6,π->4.0,δ-> .01};

(Subscript[ρ, 0]) = (nm1[0] + nm2[0])/(nm1[0] + nm2[0] + S);
(Subscript[ρ, 1]) = (nm1[0] + nm2[0])/(nm1[0] + nm2[0] + S/π);
total[t]= (nf1[t]+nf2[t]+nm1[t] + nm2[t]);
totalm[t] = (nm1[t] + nm2[t]);
nm1[t] = nf1[t];
nm2[t] = nf2[t];
Sols=NDSolve[
{
nf1'[t]== nf1[t](1-total[t]/K)(Log[b (Subscript[ρ, 0])] ) ((nm1[t]+nm2[t]/2)/totalm[t])+
nf2[t](1-total[t]/K)(Log[b (Subscript[ρ, 1])]  )(nm1[t]/(2 totalm[t])) - δ nf1[t]/.Pars,

nf2'[t]== nf2[t](1-total[t]/K)(Log[b (Subscript[ρ, 1])]  ) (((nm1[t]/2)+nm2[t])/totalm[t]) + nf1[t](1-total[t]/K)(Log[b (Subscript[ρ, 0])] )  (nm2[t]/(2 totalm[t]))- δ nf2[t]/.Pars,

freq[t]==nf2[t]/(nf1[t]+nf2[t])/.Pars,
males[t] == nf1[t] + nf2[t]/.Pars,

(*nf1[0]== 32*^4,nf2[0]== 1*^4, nm1[0] == 32*^4, nm2[0] == 1*^4*)
nf1[0]== 32.7*^4,nf2[0]== 0.3*^4, nm1[0] == 32.7*^4, nm2[0] == 0.3*^4
(*nf1[0]== 33*^4,nf2[0]== 0*^4, nm1[0] == 33*^4, nm2[0] == 0*^4*)
(*nf1[0]== 0*^4,nf2[0]== 33*^4, nm1[0] == 0*^4, nm2[0] == 33*^4*)
},

{nf1[t],nf2[t],males[t], freq[t],(Subscript[ρ, 0]),(Subscript[ρ, 1])},{t,0,1000}   ];

(*exp(r) = b (n/(n+S)), a strictly increasing function of n*)

LogPlot[{Evaluate [nf1[t]/. Sols], Evaluate[nf2[t]/.Sols](*, Evaluate[males[t]/.Sols]*) },{t,0,1000 (*1000*)},PlotRange->{1,1*^7},Frame->True,PlotStyle->{Red,{Blue ,Dashed}},PlotLegends->{"Nf1","Nf2","Males"}, FrameLabel->{Style["Time",16],Style["Density",16]}, LabelStyle->{},Epilog->{
Text[Style["non-discriminator",   24,Red,  FontFamily->Arial],{500,7.9}(*{35.3,10.3}*)] ,Text[Style["discriminator",24,Blue,FontFamily->Arial],{550,14.6}(*{65.1,4.9}*)],
Text[Style["A",24,Black,FontFamily->Arial],{950,14.5}(*{95.1,14.5}*)]
}
]
data0 = Table[{t, (Subscript[ρ, 0]),(Subscript[ρ, 1]),nf2[t]}/.Sols /.Pars,{t,1,50,24}]

Plot[{Evaluate[freq[t]/.Sols]},{t,0,(*16*)1000},PlotRange->{0,1},Frame->True,PlotStyle->{Green},PlotLegends->{"freq"},FrameLabel->{Style["Time",16],Style["Frequency",16]},Epilog->{
Text[Style["Discrimination frequency",24,FontFamily->Arial],{(*7.5*)550,0.5}],
Text[Style["B",24,FontFamily->Arial],{950,0.9}(*{15.3,0.9}*)]
}]