{
kminC=0;
kmaxC=20;
kminI=0;
kmaxI=20;

lambdaC=3.0;
lambdaI=3.0;

cc=0;
Do[{
PkC[i]=Exp[-lambdaC]*lambdaC^i/i!;
(*PkC[i]=i^-lambdaC;*)
cc=cc+PkC[i];
},{i,kminC,kmaxC}];
Do[{
PkC[i]=PkC[i]/cc;
},{i,kminC,kmaxC}];

cc=0;
Do[{
PkI[i]=Exp[-lambdaI]*lambdaI^i/i!;
cc=cc+PkI[i];
},{i,kminI,kmaxI}];
Do[{
PkI[i]=PkI[i]/cc;
},{i,kminI,kmaxI}];


kmedC=0;
Do[{
kmedC=kmedC+i*PkC[i];
},{i,kminC,kmaxC}];
kmedI=0;
Do[{
kmedI=kmedI+i*PkI[i];
},{i,kminI,kmaxI}];


G1C[x_]:=Sum[i*PkC[i]/kmedC*x^(i-1),{i,kminC,kmaxC}];
G1I[x_]:=Sum[i*PkI[i]/kmedI*x^(i-1),{i,kminI,kmaxI}];
G1T[x_]:=G1I[G1C[x]];

eps1[beta_,ff_]:=(D[G1C[beta*(1-ff)+(1-beta)*x],x]/.x->1)-(D[G1C[beta*(1-beta)*(1-ff)+(1-beta)*x],x]/.x->1);
eps2[beta_,ff_]:=(D[G1C[(1-beta)+beta*(1-ff)*x],x]/.x->1)*(D[G1T[(1-beta)+beta*x],x]/.x->1);

G1Cmed=D[G1C[x],x]/.x->1;
ListExp={};
Do[{
beta=1-0.001*i;
fi=0.5;
sol=FindRoot[(eps1[beta,ff]+eps2[beta,ff])/(beta*G1Cmed)==1,{ff,fi}];
ffsol=Abs[ff]/.sol[[1]];
If[ffsol<0.01,Break[]];(*Exit the loop if the solution becomes unstable (close to ffsol=0) *)
If[Mod[i,10]==0,{Print[i]}];
fi=ffsol;
AppendTo[ListExp,{beta,ffsol}];
},{i,1,1000}]
Export["CriticalCurve_StrategyA.dat",ListExp];

}