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

Pars={S -> 3.0*^6,b-> 10,K-> 1.0*^6,π->2.0 (*2.0 or 1.0*),m-> 1000 (*m-> 1000 or m-> 3400*),δ->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]+m/.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]+0.03125m/.Pars,

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

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

},

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

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

data0 = Table[{t,nf1[t],nf2[t],2*(nf1[t]+nf2[t])/K,freq[t]}/.Sols /.Pars,{t,100,1000,100}]

LogPlot[{Evaluate [nf1[t]/. Sols], Evaluate[nf2[t]/.Sols](*, Evaluate[males[t]/.Sols]*) },{t,0,20},PlotRange->{1,(*5*)5*^6},Frame->True,PlotStyle->{Red,{Blue ,Dashed}},PlotLegends->{"F","(SuperStar[F])"}, FrameLabel->{Style["Time",16],Style["Density",16]}, LabelStyle->{},Epilog->{
(*Inset[Framed[Style["B",16],Background->White,{19.2,12.5}]],*)
Text[Style["m = 3506" (*3506 or 1031*),24,FontFamily->Arial],{10.8,4.5}(*{20.8,4.5}*)],
Text[Style["non-discriminator",   24,Red,  FontFamily->Arial],(*{15.3,13.9}*){7.3,13.9}] ,Text[Style["discriminator",24,Blue,FontFamily->Arial],(*{13,9}*){5.1,10.0}],
Text[Style["A",24,Black,FontFamily->Arial,Background->White],{19.3,14.3}]
}]


Plot[{Evaluate[freq[t]/.Sols]},{t,0,20},PlotRange->{0,1.0},Frame->True,PlotStyle->{Green},PlotLegends->{"freq"},FrameLabel->{Style["Time",16],Style["Frequency",16]},
LabelStyle->{},Epilog->{Text[Style["Discrimination frequency",24,FontFamily->Arial],{9.5,0.5}],
Text[Style["B",24,Black,FontFamily->Arial,Background->White],{19.3,0.9}(*{39.3,0.9}*)]}]