Numerical Solutions of SIT in continuous time; here, evolution of a preference; haploid model where nmi = nfi

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

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[
{
s'[t] == 0/.Pars,
nf1'[t]== nf1[t](1-total[t]/K)(Log[b] + Log[totalm[t]] - Log[totalm[t]+s[t] ]) ((nm1[t]+nm2[t]/2)/totalm[t])+
nf2[t](1-total[t]/K)(Log[b] +  Log[π] +Log[totalm[t]] - Log[π(totalm[t])+s[t]] )(nm1[t]/(2 totalm[t]))-δ nf1[t]/.Pars,

nf2'[t]== nf2[t](1-total[t]/K)(Log[b] +  Log[π] +Log[totalm[t]] - Log[π(totalm[t])+s[t]] ) ((nm1[t]/2+nm2[t])/totalm[t]) + nf1[t](1-total[t]/K)(Log[b] + Log[totalm[t]] - Log[totalm[t]+s[t]] )  (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.7*^4,nf2[0]== .3*^4, nm1[0] == 32.7*^4, nm2[0] == .3*^4,s[0] == 3*^6

(*nf1[0]== 1.0*^5,nf2[0]== 1.0*^3, nm1[0] == 1.0*^5, nm2[0] == 1.0*^3,s[0] == 3*^6*),
(*WhenEvent[nf2[t] + nf1[t]≤ 10 && t<35, s[t] -> 00],*)
WhenEvent[nf2[t] + nf1[t]≤ 10 , s[t] -> 00],
WhenEvent[nf2[t] + nf1[t] ≥ 1000 && s[t] ≤ 200,s[t] -> 3*^6]

},

{nf1[t],nf2[t],males[t], freq[t],s[t]},{t,0,50}   ];

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

LogPlot[{Evaluate [nf1[t]/. Sols], Evaluate[nf2[t]/.Sols],Evaluate[s[t]]/.Sols},{t,0,35},PlotRange->{1,3.5*^6},Frame->True,PlotStyle->{Red,{Blue,Dashed},{Black,Thick}},FrameLabel->{Style["Time",16],Style["Density",16]}, LabelStyle->{},PlotLegends->{"F","F*","S"},Epilog->{
Text[Style["non-discriminator",   24,Red,  FontFamily->Arial],{11.4,14.0}] ,Text[Style["discriminator",24,Blue,FontFamily->Arial],{27.1,7.8}]}]

Plot[{Evaluate[freq[t]/.Sols]},{t,0,30},PlotRange->{0,1},Frame->True,PlotStyle->{Green},PlotLegends->{"freq"},FrameLabel->{Style["Time",16],Style["Frequency",16]}]