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[
{
nf1'[t]== nf1[t](1-total[t]/K)(Log[b] + Log[totalm[t]] - Log[totalm[t]+S] ) ((nm1[t]+nm2[t]/2)/totalm[t])+
nf2[t](1-total[t]/K)(Log[b] +  Log[π] +Log[totalm[t]] - Log[π(totalm[t])+S] )(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] ) ((nm1[t]/2+nm2[t])/totalm[t]) + nf1[t](1-total[t]/K)(Log[b] + Log[totalm[t]] - Log[totalm[t]+S] )  (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[t],nf2[t],males[t], freq[t]},{t,0,20}   ];

(*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,20},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],{7.3,14.3}] ,Text[Style["discriminator",24,Blue,FontFamily->Arial],{5.1,6.9}],
Text[Style["A",24,Black,FontFamily->Arial],{19.1,14.5}]
}
]

Plot[{Evaluate[freq[t]/.Sols]},{t,0,16},PlotRange->{0,1},Frame->True,PlotStyle->{Green},PlotLegends->{"freq"},FrameLabel->{Style["Time",16],Style["Frequency",16]},Epilog->{
Text[Style["Discrimination frequency",24,FontFamily->Arial],{6.5,0.5}],
Text[Style["B",24,FontFamily->Arial],{15.1,0.9}]
}]