Calculations for Two Nested Limit Cycles in Two-Species Reactions with Wolfram Mathematica

Created by Ilona Nagy, Valery Romanovski and János Tóth
Last modified: 24 August, 2020

Model 1  with a stable outer and an unstable inner limit cycle

Preparations

System (3)

x'=(x^2) y+x y-c1 (x^2)-d1 x+e1 y+f1
y'=-(x^2)y-x y+c1 (x^2)+d2 x-e2 y+f2

c1,d1,d2,e1,e2,f1,f2≥0

The ReactionKinetics program package

The ReactionKinetics program package is available at http://extras.springer.com (ISBN: 978-1-4939-8643-9).
It can be used if either ReactionKinetics.m is put in the same folder as this notebook or ReactionKinetics.wl is added in the packages in the applications.

Quit

SetDirectory[NotebookDirectory[]];
SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{ContourPlot,DateListPlot,Plot,ListLinePlot,ListPlot,ListLogPlot,LogLinearPlot,LogPlot,ParametricPlot,Plot3D,RegionPlot};
LaunchKernels[];
Needs["ReactionKinetics`"];

Figure 1: Creating the reaction graph with the help of the ReactionKinetics package

ClearAll[model1];model1={2X+Y->3X,X->0,2X->X+Y->2X,0->X->X+Y,0->Y->X+Y,Y-> 0};
RightHandSide[{model1},{1,d1,c1,1,f1,d2,f2,e1,e2},{x,y}]
model1fig=ShowFHJGraph[model1,{1,(Subscript[d, 1]),(Subscript[c, 1]),1,(Subscript[f, 1]),(Subscript[d, 2]),(Subscript[f, 2]),(Subscript[e, 1]),(Subscript[e, 2])},DirectedEdges->True, 
 VertexLabels->Automatic, EdgeLabelStyle->Directive[Red,Italic,16],VertexLabelStyle->Directive[Blue,16],GraphLayout->"TutteEmbedding"]

Export["Fig-1-Model1.pdf",model1fig]

The singular point is shifted into the origin

Equations (4)-(5): Singular points if x0=1

Quit

ClearAll[xd,yd,x,y,c1,d1,d2,e1,e2,f1,f2];
xd=(x^2) y+x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y-x y+c1 (x^2)+d2 x-e2 y+f2;
Solve[{xd==0,yd==0}/.x->1,{d1,y}]//FullSimplify

Equation (6): The singular point (if x0=1) is shifted into (0,0)

ClearAll[xd,yd,x,y,c1,d1,d2,e1,e2,f1,f2,x0,y0,x1d,y1d,xx1,yy1];
xd=(x^2) y+x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y-x y+c1 (x^2)+d2 x-e2 y+f2;
x0=1;
d1=((d2 (2+e1)+c1 (e1-e2)+(2+e2) f1+(2+e1) f2)/(2+e2));y0=((c1+d2+f2)/(2+e2));
xx1=x-x0;yy1=y-y0;
x1d=D[xx1,x]xd+D[xx1,y]yd/.{x->x1+x0,y->y1+y0}//Factor
y1d=D[yy1,x]xd+D[yy1,y]yd/.{x->x1+x0,y->y1+y0}//Factor

The Jacobian at the origin

Equations (7)-(8)

Jac=D[{x1d,y1d},{{x1,y1}}];
JacOrigin=Jac/.{x1->0,y1->0}//Simplify

trace=Tr[JacOrigin]//Factor
Solve[trace==0,c1]//FullSimplify

c1=((d2-d2 e1-(2+e2) (2+e2+f1)+f2-e1 f2)/(1+e1+e2));
evalues=Eigenvalues[JacOrigin]//FullSimplify

Equation (9)

beta=-((evalues[[1]])^2)//Factor

System (10)

pp=x1d/.{x1->x,y1->y}//Factor
qq=y1d/.{x1->x,y1->y}//Factor

Lyapunov’s theorem

System (10)

Quit

pp=(1/(1+e1+e2))(2 x+2 e1 x+3 e2 x+e1 e2 x+(e2^2) x+2 (x^2)+d2 e1 (x^2)+3 e2 (x^2)+(e2^2) (x^2)+f1 (x^2)+e2 f1 (x^2)+e1 f2 (x^2)+2 y+3 e1 y+(e1^2) y+2 e2 y+e1 e2 y+3 x y+3 e1 x y+3 e2 x y+(x^2) y+e1 (x^2) y+e2 (x^2) y);
qq=-(1/(1+e1+e2))(2 x+d2 e1 x+5 e2 x-d2 e2 x+2 (e2^2) x+f1 x+2 e2 f1 x+f2 x+2 e1 f2 x+2 (x^2)+d2 e1 (x^2)+3 e2 (x^2)+(e2^2) (x^2)+f1 (x^2)+e2 f1 (x^2)+e1 f2 (x^2)+2 y+2 e1 y+3 e2 y+e1 e2 y+(e2^2) y+3 x y+3 e1 x y+3 e2 x y+(x^2) y+e1 (x^2) y+e2 (x^2) y);

Program

Ser[s_]:=Plus@@Table[(x^i)(y^(s-i))p[i,s-i],{i,0,s}];
Hom[s_]:=Table[p[s-i,i],{i,0,s}];
hh=Sum[Ser[i],{i,2,6}];(*9*)
Lie=D[hh,x] pp+D[hh,y] qq//Expand;
RHS=g1((((x^2)+(y^2)))^2)+g2((((x^2)+(y^2)))^3)+g3((((x^2)+(y^2)))^4)//Expand;
vv=Lie-RHS//Expand;
CoefPol[f_,s_]:=Module[{m,lis,t},lis={}; m=Expand[f];  Do[Do[If[i+j==s,lis=AppendTo[lis,Coefficient[m,x^i y^j ]/.{x->0,y->0,z->0}]],{i,0,s}],{j,0,s}];
ls[s]=lis];
Do[CoefPol[vv,i],{i,1,9}];

Degree 1, 2

ls[1]
ls[2]//Factor;
sol2=Solve[ls[2]==0,Hom[2]]//Simplify;
{p[2,0],p[1,1],p[0,2]}={p[2,0],p[1,1],p[0,2]}/.sol2[[1]];
ls[2]//Simplify

Quadratic form

ClearAll[qv,mat,a11,det,eg];
qv=Ser[2]//FullSimplify
mat=1/2 D[qv,{{x,y},2}];
mat//MatrixForm
a11=mat[[1,1]]
det=Det[mat]//Factor
eg=Eigenvalues[mat]//Simplify//Factor;

Equation (13)

p[1,1]=1;
qv
mat//MatrixForm
a11
det

Conditions for a positive definite quadratic form

d1=((d2 (2+e1)+c1 (e1-e2)+(2+e2) f1+(2+e1) f2)/(2+e2));y0=((c1+d2+f2)/(2+e2));
c1=((d2-d2 e1-(2+e2) (2+e2+f1)+f2-e1 f2)/(1+e1+e2));
beta=(1/(1+e1+e2))(-2 e1+2 d2 e1+d2 (e1^2)+2 e2-2 d2 e2+e1 e2-d2 e1 e2-(e2^2)+e1 (e2^2)-(e2^3)+2 f1+e1 f1+4 e2 f1+2 e1 e2 f1+2 f2+5 e1 f2+2 (e1^2) f2);

Reduce[ a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0&&f1>0 &&f2>0,{d2,e1,e2,f1,f2} ]//FullSimplify

Setting f1 and f2: Equation (14)

ClearAll[f1,f2];
f1=1;f2=2;

Reduce[ a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e1,e2} ]//FullSimplify

Degree 3

ls[3]//Factor;
sol3=Solve[ls[3]==0,Hom[3]]//Simplify;
{p[3,0],p[2,1],p[1,2],p[0,3]}={p[3,0],p[2,1],p[1,2],p[0,3]}/.sol3[[1]]//Factor;
ls[3]//Simplify

Degree 4: Equation (15)

ls[4]//Factor;
sol4=Solve[ls[4]==0,AppendTo[Hom[4],g1]];
{g1,p[4,0],p[3,1],p[2,2],p[1,3],p[0,4]}={g1,p[4,0],p[3,1],p[2,2],p[1,3],p[0,4]}/.sol4[[1]];
g1//Factor

Variables[g1]

Degree 5

ls[5]//Factor;
sol5=Solve[ls[5]==0,Hom[5]];
{p[5,0],p[4,1],p[3,2],p[2,3],p[1,4],p[0,5]}={p[5,0],p[4,1],p[3,2],p[2,3],p[1,4],p[0,5]}/.sol5[[1]]//Factor;
ls[5]//Factor

Degree 6

ls[6]//Factor;
sol6=Solve[ls[6]==0,AppendTo[Hom[6],g2]];
{p[6,0],p[5,1],p[4,2],p[3,3],p[2,4],p[1,5],p[0,6],g2}={p[6,0],p[5,1],p[4,2],p[3,3],p[2,4],p[1,5],p[0,6],g2}/.sol6[[1]]//Factor;

g2
Variables[g2]

p[2,2]=0;
g2

Variables[g2]

Finding possible values of e1 such that g1=0

ClearAll[d2,e1,e2];e1=(1/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify

ClearAll[d2,e1,e2];e1=(2/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify

ClearAll[d2,e1,e2];e1=(3/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify

ClearAll[d2,e1,e2];e1=(4/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify

ClearAll[d2,e1,e2];e1=(5/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify

ClearAll[d2,e1,e2];e1=(6/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify

ClearAll[d2,e1,e2];e1=(7/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify

ClearAll[d2,e1,e2];e1=(8/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify

ClearAll[d2,e1,e2];e1=(9/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify

ClearAll[d2,e1,e2];e1=1;
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e2} ]//FullSimplify

Setting the value of e1, d2: Equations (16)-(18)

ClearAll[d2,e1,e2];e1=(5/10);d2=20;
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&e2>0,{e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&e2>0,{e2} ]//FullSimplify

g1//Factor

g2//Factor

e2=(NumericalApproximation[Root, (RowDefault["0.771", FEExpressionsNumericalApproximationElider]), Root[-14027+15420 #1+2826 ((#1)^2)+926 ((#1)^3)+72 ((#1)^4)+8 ((#1)^5)&,1], ]);
e2//N

Perturbation

Setting e2 such that g1=0,g2<0, trace(J)=0 (no perturbation)

Quit

ClearAll[xd,yd,c1,d1,d2,e1,e2,f1,f2,g1,g2,sol];
xd=(x^2) y+x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y-x y+c1 (x^2)+d2 x-e2 y+f2;
f1=1;f2=2;
e1=(5/10);d2=20;
d1=((d2 (2+e1)+c1 (e1-e2)+(2+e2) f1+(2+e1) f2)/(2+e2));
c1=((d2-d2 e1-(2+e2) (2+e2+f1)+f2-e1 f2)/(1+e1+e2));
e2=(NumericalApproximation[Root, (RowDefault["0.771", FEExpressionsNumericalApproximationElider]), Root[-14027+15420 #1+2826 ((#1)^2)+926 ((#1)^3)+72 ((#1)^4)+8 ((#1)^5)&,1], ]);
g1=((2 (-14027+15420 e2+2826 (e2^2)+926 (e2^3)+72 (e2^4)+8 (e2^5)))/((2+e2) (17163-19172 e2+12044 (e2^2)-1888 (e2^3)+256 (e2^4))));
g2=-((2 (-126512084933352295+478399692835658985 e2-762458375238510816 (e2^2)+730329622432012315 (e2^3)-466875108129002500 (e2^4)+180057388265535584 (e2^5)-42412278548678749 (e2^6)+7565318705371668 (e2^7)-1899913140260206 (e2^8)+200501244715092 (e2^9)+44607817747872 (e2^10)+18408318367872 (e2^11)+8352995841088 (e2^12)+1190085589888 (e2^13)+348608893184 (e2^14)+33760183296 (e2^15)+6787676672 (e2^16)+473953280 (e2^17)+64299008 (e2^18)+2555904 (e2^19)+196608 (e2^20)))/(3 (2+e2) (83-42 e2+8 (e2^2)) (((-73+85 e2+(e2^2)+2 (e2^3)))^2) (109-25 e2+32 (e2^2)+4 (e2^3)) (17163-19172 e2+12044 (e2^2)-1888 (e2^3)+256 (e2^4)) (23933-29628 e2+23764 (e2^2)-2976 (e2^3)+512 (e2^4))));
{g1,g2}//N

ClearAll[sol];
sol=Solve[{xd==0,yd==0},{x,y},Reals][[1]]//N
D[{xd,yd},{{x,y}}]/.sol;
Eigenvalues[D[{xd,yd},{{x,y}}]/.sol]

3.2.1. Perturbing e2 such that g1>0,g2<0, trace(J)=0

Quit

ClearAll[xd,yd,c1,d1,d2,e1,e2,f1,f2,g1,g2,sol];
xd=(x^2) y+x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y-x y+c1 (x^2)+d2 x-e2 y+f2;
f1=1;f2=2;
e1=(5/10);d2=20;
d1=((d2 (2+e1)+c1 (e1-e2)+(2+e2) f1+(2+e1) f2)/(2+e2));
c1=((d2-d2 e1-(2+e2) (2+e2+f1)+f2-e1 f2)/(1+e1+e2));
e2=(78/100);(*perturbed parameter for g1>0*)
g1=((2 (-14027+15420 e2+2826 (e2^2)+926 (e2^3)+72 (e2^4)+8 (e2^5)))/((2+e2) (17163-19172 e2+12044 (e2^2)-1888 (e2^3)+256 (e2^4))));
g2=-((2 (-126512084933352295+478399692835658985 e2-762458375238510816 (e2^2)+730329622432012315 (e2^3)-466875108129002500 (e2^4)+180057388265535584 (e2^5)-42412278548678749 (e2^6)+7565318705371668 (e2^7)-1899913140260206 (e2^8)+200501244715092 (e2^9)+44607817747872 (e2^10)+18408318367872 (e2^11)+8352995841088 (e2^12)+1190085589888 (e2^13)+348608893184 (e2^14)+33760183296 (e2^15)+6787676672 (e2^16)+473953280 (e2^17)+64299008 (e2^18)+2555904 (e2^19)+196608 (e2^20)))/(3 (2+e2) (83-42 e2+8 (e2^2)) (((-73+85 e2+(e2^2)+2 (e2^3)))^2) (109-25 e2+32 (e2^2)+4 (e2^3)) (17163-19172 e2+12044 (e2^2)-1888 (e2^3)+256 (e2^4)) (23933-29628 e2+23764 (e2^2)-2976 (e2^3)+512 (e2^4))));
{c1,d1}
{c1,d1}//N
{g1,g2}//N

ClearAll[sol];
sol=Solve[{xd==0,yd==0},{x,y},Reals][[1]]//N
D[{xd,yd},{{x,y}}]/.sol;
Eigenvalues[D[{xd,yd},{{x,y}}]/.sol]//FullSimplify

3.2.2. Perturbing c1 and e2 such that g1>0,g2<0, trace(J)<0

Quit

ClearAll[xd,yd,c1,d1,d2,e1,e2,f1,f2,g1,g2,sol];
xd=(x^2) y+x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y-x y+c1 (x^2)+d2 x-e2 y+f2;
f1=1;f2=2;
e1=(5/10);d2=20;
d1=((d2 (2+e1)+c1 (e1-e2)+(2+e2) f1+(2+e1) f2)/(2+e2));
c1=(22/100);(*perturbed parameter for trace(J)<0*)
e2=(78/100);(*perturbed parameter for g1>0*)
g1=((2 (-14027+15420 e2+2826 (e2^2)+926 (e2^3)+72 (e2^4)+8 (e2^5)))/((2+e2) (17163-19172 e2+12044 (e2^2)-1888 (e2^3)+256 (e2^4))));
g2=-((2 (-126512084933352295+478399692835658985 e2-762458375238510816 (e2^2)+730329622432012315 (e2^3)-466875108129002500 (e2^4)+180057388265535584 (e2^5)-42412278548678749 (e2^6)+7565318705371668 (e2^7)-1899913140260206 (e2^8)+200501244715092 (e2^9)+44607817747872 (e2^10)+18408318367872 (e2^11)+8352995841088 (e2^12)+1190085589888 (e2^13)+348608893184 (e2^14)+33760183296 (e2^15)+6787676672 (e2^16)+473953280 (e2^17)+64299008 (e2^18)+2555904 (e2^19)+196608 (e2^20)))/(3 (2+e2) (83-42 e2+8 (e2^2)) (((-73+85 e2+(e2^2)+2 (e2^3)))^2) (109-25 e2+32 (e2^2)+4 (e2^3)) (17163-19172 e2+12044 (e2^2)-1888 (e2^3)+256 (e2^4)) (23933-29628 e2+23764 (e2^2)-2976 (e2^3)+512 (e2^4))));
{c1,d1}
{c1,d1}//N
{g1,g2}//N

ClearAll[sol];
sol=Solve[{xd==0,yd==0},{x,y},Reals][[1]]//N
D[{xd,yd},{{x,y}}]/.sol;
Eigenvalues[D[{xd,yd},{{x,y}}]/.sol]

-0.0017985611510798125`*2

3.2.1. Plotting the limit cycles when  g1>0,g2<0, trace(J)=0

Preparations

Quit

SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ParametricPlot,ListPlot,ListLinePlot};
SetDirectory[NotebookDirectory[]];
SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ListPlot,ListLinePlot,ListLogLogPlot,ParametricPlot,DateListPlot,DiscretePlot};
LaunchKernels[];

The function creating the plots

ClearAll[p,q,x,y,c1,d1,d2,e1,e2,f1,f2];
p[x_,y_]:=(x^2) y+x y-c1 (x^2)-d1 x+e1 y+f1;
q[x_,y_]:=-(x^2)y-x y+c1 (x^2)+d2 x-e2 y+f2;
f1=1;f2=2;
e1=(5/10);d2=20;
d1=((d2 (2+e1)+c1 (e1-e2)+(2+e2) f1+(2+e1) f2)/(2+e2));
c1=((d2-d2 e1-(2+e2) (2+e2+f1)+f2-e1 f2)/(1+e1+e2));
e2=(78/100);(*perturbed parameter for g1>0*)
{e2,e1,d2,d1,f1,f2,c1}
{e2,e1,d2,d1,f1,f2,c1}//N

ClearAll[nsol,ev,plotter];nsol=First@NSolve[Join@@Thread/@{{p[x,y],q[x,y]}==0,{x,y}>0},{x,y},20]
ev=Eigenvalues[D[{p[x,y],q[x,y]},{{x,y}}]/.nsol]
plotter[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:1000,ar_:Automatic,opts___]:=Module[{startingpoint,sys,solution,plot1},
startingpoint=({x,y}/.nsol)+shift;
sys:=NDSolveValue[Join[{u'[t]==p[u[t],v[t]],v'[t]==q[u[t],v[t]]},Thread[{u[0],v[0]}==startingpoint]],
{u,v},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
solution[t_]:=Through[sys[t]];
{ParametricPlot[Evaluate[solution[t]],{t,0,τ},Epilog->{Red,PointSize[0.05],Point[startingpoint],
Orange,Point[{x,y}/.nsol]},PlotRange->All,PlotPoints->pp,AspectRatio->ar,AxesLabel->{x,y},LabelStyle->Directive[14],ImageSize->200],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,x},LabelStyle->Directive[12],ImageSize->200],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,y},LabelStyle->Directive[12],ImageSize->200]}]

Plotter with arrow

ClearAll[nsol,ev,plotterarrow];nsol=First@NSolve[Join@@Thread/@{{p[x,y],q[x,y]}==0,{x,y}>0},{x,y},20]
ev=Eigenvalues[D[{p[x,y],q[x,y]},{{x,y}}]/.nsol]
plotterarrow[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:1000,ar_:Automatic,arrow_,opts___]:=Module[{startingpoint,sys,solution,plot1},
startingpoint=({x,y}/.nsol)+shift;
sys:=NDSolveValue[Join[{u'[t]==p[u[t],v[t]],v'[t]==q[u[t],v[t]]},Thread[{u[0],v[0]}==startingpoint]],
{u,v},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
solution[t_]:=Through[sys[t]];
{ParametricPlot[Evaluate[solution[t]],{t,0,τ},Epilog->{Black,Arrowheads->0.07,Arrow[{startingpoint,arrow}],Red,PointSize[0.05],Point[startingpoint],Orange,
 Point[{x,y}/.nsol]
},PlotRange->All,PlotPoints->pp,AspectRatio->ar,AxesLabel->{x,y},LabelStyle->Directive[14],ImageSize->200],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,x},LabelStyle->Directive[12],ImageSize->200],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,y},LabelStyle->Directive[12],ImageSize->200]}]

Figure 2

plotter[100,{0,10},Automatic,100,1000,Automatic,Method->"BDF"]
plotterarrow[100,{0,10},Automatic,100,1000,Automatic,{3.586,15.36},Method->"BDF"]

Figure 3

plotter[120,{0,0.038},Automatic,100,1000,Automatic,Method->"BDF"]

3.2.2. Plotting the limit cycles when  g1>0,g2<0, trace(J)<0

Preparations

Quit

SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ParametricPlot,ListPlot,ListLinePlot};
SetDirectory[NotebookDirectory[]];
SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ListPlot,ListLinePlot,ListLogLogPlot,ParametricPlot,DateListPlot,DiscretePlot};
LaunchKernels[];

The function creating the plots

ClearAll[p,q,x,y,c1,d1,d2,e1,e2,f1,f2];
p[x_,y_]:=(x^2) y+x y-c1 (x^2)-d1 x+e1 y+f1;
q[x_,y_]:=-(x^2)y-x y+c1 (x^2)+d2 x-e2 y+f2;
f1=1;f2=2;
e1=(5/10);d2=20;
d1=((d2 (2+e1)+c1 (e1-e2)+(2+e2) f1+(2+e1) f2)/(2+e2));
c1=(22/100);(*perturbed parameter for trace(J)<0*)
e2=(78/100);(*perturbed parameter for g1>0*)
{e2,e1,d2,d1,f1,f2,c1}
{e2,e1,d2,d1,f1,f2,c1}//N

ClearAll[nsol,ev,plotter];nsol=First@NSolve[Join@@Thread/@{{p[x,y],q[x,y]}==0,{x,y}>0},{x,y},20]
ev=Eigenvalues[D[{p[x,y],q[x,y]},{{x,y}}]/.nsol]
plotter[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:1000,ar_:Automatic,opts___]:=Module[{startingpoint,sys,solution,plot1},
startingpoint=({x,y}/.nsol)+shift;
sys:=NDSolveValue[Join[{u'[t]==p[u[t],v[t]],v'[t]==q[u[t],v[t]]},Thread[{u[0],v[0]}==startingpoint]],
{u,v},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
solution[t_]:=Through[sys[t]];
{ParametricPlot[Evaluate[solution[t]],{t,0,τ},Epilog->{Red,PointSize[0.05],Point[startingpoint],Orange, Point[{x,y}/.nsol]},PlotRange->All,PlotPoints->pp,AspectRatio->ar,AxesLabel->{x,y},LabelStyle->Directive[14],ImageSize->200],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,x},LabelStyle->Directive[12],ImageSize->200],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,y},LabelStyle->Directive[12],ImageSize->200]}]

Plotter with arrow

ClearAll[nsol,ev,plotterarrow];nsol=First@NSolve[Join@@Thread/@{{p[x,y],q[x,y]}==0,{x,y}>0},{x,y},20]
ev=Eigenvalues[D[{p[x,y],q[x,y]},{{x,y}}]/.nsol]
plotterarrow[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:1000,ar_:Automatic,arrow_,opts___]:=Module[{startingpoint,sys,solution,plot1},
startingpoint=({x,y}/.nsol)+shift;
sys:=NDSolveValue[Join[{u'[t]==p[u[t],v[t]],v'[t]==q[u[t],v[t]]},Thread[{u[0],v[0]}==startingpoint]],
{u,v},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
solution[t_]:=Through[sys[t]];
{ParametricPlot[Evaluate[solution[t]],{t,0,τ},Epilog->{Black,Arrowheads->0.07,Arrow[{startingpoint,arrow}],Red,PointSize[0.05],Point[startingpoint],Orange,
 Point[{x,y}/.nsol]
},PlotRange->All,PlotPoints->pp,AspectRatio->ar,AxesLabel->{x,y},LabelStyle->Directive[14],ImageSize->200],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,x},LabelStyle->Directive[12],ImageSize->200],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,y},LabelStyle->Directive[12],ImageSize->200]}]

Figure 4

plotter[100,{0,10},Automatic,100,1000,Automatic,Method->"BDF"]
plotterarrow[100,{0,10},Automatic,100,1000,Automatic,{3.586,15.36},Method->"BDF"]

Figure 5

plotter[120,{0,0.042},Automatic,100,1000,Automatic,Method->"BDF"]

Figure 6

plotter[120,{0,0.038},Automatic,100,1000,Automatic,Method->"BDF"]
plotterarrow[120,{0,0.038},Automatic,100,1000,Automatic,{1.066,7.956},Method->"BDF"]

Figure 7

plotter[200,{0,0.01},Automatic,100,1000,Automatic,Method->"BDF"]
plotterarrow[200,{0,0.01},Automatic,100,1000,Automatic,{1.007,7.994},Method->"BDF"]

Model 2 with an unstable outer a stable inner limit cycle

Preparations

System (19)

x'=(x^2) y-x y-c1 (x^2)-d1 x+e1 y+f1
y'=-(x^2)y+x y+c1 (x^2)+d2 x-e2 y+f2

c1,d1,d2,e1,e2,f1,f2≥0

The ReactionKinetics program package

The ReactionKinetics program package is available at http://extras.springer.com (ISBN: 978-1-4939-8643-9).
It can be used if either ReactionKinetics.m is put in the same folder as this notebook or ReactionKinetics.wl is added in the packages in the applications.

Quit

SetDirectory[NotebookDirectory[]];
SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{ContourPlot,DateListPlot,Plot,ListLinePlot,ListPlot,ListLogPlot,LogLinearPlot,LogPlot,ParametricPlot,Plot3D,RegionPlot};
LaunchKernels[];
Needs["ReactionKinetics`"];

Figure 8: Creating the reaction graph with the help of the ReactionKinetics package

ClearAll[model2];model2={2X+Y->3X,X->0,2X->X+Y->2Y,0->X->X+Y,0->Y->X+Y,Y-> 0};
RightHandSide[{model2},{1,d1,c1,1,f1,d2,f2,e1,e2},{x,y}]
model2fig=ShowFHJGraph[model2,{1,(Subscript[d, 1]),(Subscript[c, 1]),1,(Subscript[f, 1]),(Subscript[d, 2]),(Subscript[f, 2]),(Subscript[e, 1]),(Subscript[e, 2])},DirectedEdges->True, 
 VertexLabels->Automatic, EdgeLabelStyle->Directive[Red,Italic,16],VertexLabelStyle->Directive[Blue,16],GraphLayout->"TutteEmbedding"]

Export["Fig-8-Model2.pdf",model2fig]

The singular point is shifted into the origin

Singular points if x0=1

Quit

ClearAll[xd,yd,x,y,c1,d1,d2,e1,e2,f1,f2];
xd=(x^2) y-x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y+x y+c1 (x^2)+d2 x-e2 y+f2;
Solve[{xd==0,yd==0}/.x->1,{d1,y}]//FullSimplify

The singular point (if x0=1) is shifted into (0,0)

ClearAll[xd,yd,x,y,c1,d1,d2,e1,e2,f1,f2,x0,y0,x1d,y1d,xx1,yy1];
xd=(x^2) y-x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y+x y+c1 (x^2)+d2 x-e2 y+f2;
x0=1;
d1=((c1 (e1-e2)+e2 f1+e1 (d2+f2))/e2);y0=((c1+d2+f2)/e2);
xx1=x-x0;yy1=y-y0;
x1d=D[xx1,x]xd+D[xx1,y]yd/.{x->x1+x0,y->y1+y0}//Factor
y1d=D[yy1,x]xd+D[yy1,y]yd/.{x->x1+x0,y->y1+y0}//Factor

The Jacobian at the origin

Jac=D[{x1d,y1d},{{x1,y1}}];
JacOrigin=Jac/.{x1->0,y1->0}//Simplify

trace=Tr[JacOrigin]//Factor
Solve[trace==0,c1]//FullSimplify

c1=((d2-d2 e1-e2 (e2+f1)+f2-e1 f2)/(-1+e1+e2));
evalues=Eigenvalues[JacOrigin]//FullSimplify

beta=-((evalues[[1]])^2)//Factor

pp=x1d/.{x1->x,y1->y}//Factor
qq=y1d/.{x1->x,y1->y}//Factor

Lyapunov’s theorem

System

Quit

pp=(1/(-1+e1+e2))(-e2 x+e1 e2 x+(e2^2) x+d2 e1 (x^2)-e2 (x^2)+(e2^2) (x^2)-f1 (x^2)+e2 f1 (x^2)+e1 f2 (x^2)-e1 y+(e1^2) y+e1 e2 y-x y+e1 x y+e2 x y-(x^2) y+e1 (x^2) y+e2 (x^2) y);
qq=-(1/(-1+e1+e2))(d2 e1 x-e2 x-d2 e2 x+2 (e2^2) x-f1 x+2 e2 f1 x-f2 x+2 e1 f2 x+d2 e1 (x^2)-e2 (x^2)+(e2^2) (x^2)-f1 (x^2)+e2 f1 (x^2)+e1 f2 (x^2)-e2 y+e1 e2 y+(e2^2) y-x y+e1 x y+e2 x y-(x^2) y+e1 (x^2) y+e2 (x^2) y);

Program

Ser[s_]:=Plus@@Table[(x^i)(y^(s-i))p[i,s-i],{i,0,s}];
Hom[s_]:=Table[p[s-i,i],{i,0,s}];
hh=Sum[Ser[i],{i,2,6}];(*9*)
Lie=D[hh,x] pp+D[hh,y] qq//Expand;
RHS=g1((((x^2)+(y^2)))^2)+g2((((x^2)+(y^2)))^3)+g3((((x^2)+(y^2)))^4)//Expand;
vv=Lie-RHS//Expand;
CoefPol[f_,s_]:=Module[{m,lis,t},lis={}; m=Expand[f];  Do[Do[If[i+j==s,lis=AppendTo[lis,Coefficient[m,x^i y^j ]/.{x->0,y->0,z->0}]],{i,0,s}],{j,0,s}];
ls[s]=lis];
Do[CoefPol[vv,i],{i,1,9}];

Degree 1, 2

ls[1]
ls[2]//Factor;
sol2=Solve[ls[2]==0,Hom[2]]//Simplify;
{p[2,0],p[1,1],p[0,2]}={p[2,0],p[1,1],p[0,2]}/.sol2[[1]];
ls[2]//Simplify

Quadratic form

ClearAll[qv,mat,a11,det,eg];
qv=Ser[2]//FullSimplify
mat=1/2 D[qv,{{x,y},2}];
mat//MatrixForm
a11=mat[[1,1]]
det=Det[mat]//Factor
eg=Eigenvalues[mat]//Simplify//Factor;

Conditions for a positive definite quadratic form

d1=((c1 (e1-e2)+e2 f1+e1 (d2+f2))/e2);y0=((c1+d2+f2)/e2);
c1=((d2-d2 e1-e2 (e2+f1)+f2-e1 f2)/(-1+e1+e2));
beta=(1/(-1+e1+e2))(d2 (e1^2)-e1 e2-d2 e1 e2+(e2^2)+e1 (e2^2)-(e2^3)-e1 f1+2 e1 e2 f1-e1 f2+2 (e1^2) f2);

Reduce[ a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0&&f1>0 &&f2>0,{d2,e1,e2,f1,f2} ]//FullSimplify

Setting f1 and f2

ClearAll[f1,f2];
p[1,1]=1;f1=(1/2);f2=(1/10);

Reduce[ a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e1,e2} ]//FullSimplify

Degree 3

ls[3]//Factor;
sol3=Solve[ls[3]==0,Hom[3]]//Simplify;
{p[3,0],p[2,1],p[1,2],p[0,3]}={p[3,0],p[2,1],p[1,2],p[0,3]}/.sol3[[1]]//Factor;
ls[3]//Simplify

Degree 4

ls[4]//Factor;
sol4=Solve[ls[4]==0,AppendTo[Hom[4],g1]];
{g1,p[4,0],p[3,1],p[2,2],p[1,3],p[0,4]}={g1,p[4,0],p[3,1],p[2,2],p[1,3],p[0,4]}/.sol4[[1]];
g1//Factor

Variables[g1]

Degree 5

ls[5]//Factor;
sol5=Solve[ls[5]==0,Hom[5]];
{p[5,0],p[4,1],p[3,2],p[2,3],p[1,4],p[0,5]}={p[5,0],p[4,1],p[3,2],p[2,3],p[1,4],p[0,5]}/.sol5[[1]]//Factor;
ls[5]//Factor

Degree 6

ls[6]//Factor;
sol6=Solve[ls[6]==0,AppendTo[Hom[6],g2]];
{p[6,0],p[5,1],p[4,2],p[3,3],p[2,4],p[1,5],p[0,6],g2}={p[6,0],p[5,1],p[4,2],p[3,3],p[2,4],p[1,5],p[0,6],g2}/.sol6[[1]]//Factor;

g2
Variables[g2]

p[2,2]=0;
g2;

Variables[g2]

Investigating the possible values of e1 and d2

ClearAll[d2,e1,e2];e1=(1/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2==0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify

ClearAll[d2,e1,e2];e1=(2/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2==0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify

ClearAll[d2,e1,e2];e1=(3/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2==0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify

Setting the value of d2

ClearAll[d2,e1,e2];d2=1/100;
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&e1≥0&&e2>0,{e1,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&e1≥0&&e2>0,{e1,e2} ]//FullSimplify
Reduce[ g1==0&&g2==0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&e1≥0&&e2>0,{e1,e2} ]//FullSimplify

Setting the value of e1, d2 such that g1=0, g2>0

ClearAll[d2,e1,e2];e1=(18/100);d2=(1/100);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&e2>0,{e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&e2>0,{e2} ]//FullSimplify
Reduce[ g1==0&&g2==0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&e2>0,{e2} ]//FullSimplify

g1//Factor

g2//Factor

e2=((NumericalApproximation[Root, (RowDefault["266.", FEExpressionsNumericalApproximationElider]), Root[-107789454217728+604841338944 #1-502933212 ((#1)^2)-1266174 ((#1)^3)+988 ((#1)^4)+((#1)^5)&,3], ])/900);
e2//N

Perturbation

1. Setting e2 such that g1=0,g2>0, trace(J)=0

Quit

ClearAll[xd,yd,c1,d1,d2,e1,e2,f1,f2,g1,g2,sol];
xd=(x^2) y-x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y+x y+c1 (x^2)+d2 x-e2 y+f2;
f1=(1/2);f2=(1/10);
e1=(18/100);d2=(1/100);
d1=((c1 (e1-e2)+e2 f1+e1 (d2+f2))/e2);
c1=((d2-d2 e1-e2 (e2+f1)+f2-e1 f2)/(-1+e1+e2));
e2=((NumericalApproximation[Root, (RowDefault["266.", FEExpressionsNumericalApproximationElider]), Root[-107789454217728+604841338944 #1-502933212 ((#1)^2)-1266174 ((#1)^3)+988 ((#1)^4)+((#1)^5)&,3], ])/900);
g1=(4 (-256700232+1296384900 e2-970164375 (e2^2)-2198218750 (e2^3)+1543750000 (e2^4)+1406250000 (e2^5)))/(625 e2 (29488131-8127900 e2-113832500 (e2^2)-149000000 (e2^3)+400000000 (e2^4)));
g2=-((2 (-2134161091674039218421625601149503+140355978576685509466108656363508650 e2-1754283565687266330660162078035242500 (e2^2)+8737367532700018023712655697839625000 (e2^3)-21589750138203633706431942941437500000 (e2^4)+75104078381422305038823485272500000000 (e2^5)-513431090081409858991186790953125000000 (e2^6)+2056162232816920023362338942968750000000 (e2^7)-3838597863908140206223129570312500000000 (e2^8)+1617903538459282225849587890625000000000 (e2^9)+4615310325659104860692578125000000000000 (e2^10)-1399763372530238358330078125000000000000 (e2^11)-20086555324348486795898437500000000000000 (e2^12)+40652986288855743115234375000000000000000 (e2^13)-41911717827629003906250000000000000000000 (e2^14)+37465861770189697265625000000000000000000 (e2^15)-37513628819775390625000000000000000000000 (e2^16)+29057543024902343750000000000000000000000 (e2^17)-12113875488281250000000000000000000000000 (e2^18)+1593969726562500000000000000000000000000 (e2^19)+263671875000000000000000000000000000000 (e2^20)))/(1875 e2 (-3549+850 e2+10000 (e2^2)) (((25299+450 e2-295000 (e2^2)+250000 (e2^3)))^2) (-25299-450 e2-320000 (e2^2)+500000 (e2^3)) (29488131-8127900 e2-113832500 (e2^2)-149000000 (e2^3)+400000000 (e2^4)) (38082789-250500 e2-60467500 (e2^2)-515000000 (e2^3)+800000000 (e2^4))));
{c1,d1,e2}//N
{g1,g2}//N

Solve[{xd==0,yd==0},{x,y},Reals]//N

ClearAll[sol];
sol=Solve[{xd==0,yd==0},{x,y},Reals][[3]]//N
D[{xd,yd},{{x,y}}]/.sol;
Eigenvalues[D[{xd,yd},{{x,y}}]/.sol]//N

2. Perturbing e2 such that g1<0,g2>0, trace(J)=0

Quit

ClearAll[xd,yd,c1,d1,d2,e1,e2,f1,f2,g1,g2,sol];
xd=(x^2) y-x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y+x y+c1 (x^2)+d2 x-e2 y+f2;
f1=(1/2);f2=(1/10);
e1=(18/100);d2=(1/100);
d1=((c1 (e1-e2)+e2 f1+e1 (d2+f2))/e2);
c1=((d2-d2 e1-e2 (e2+f1)+f2-e1 f2)/(-1+e1+e2));
e2=(278/1000);(*perturbed parameter for g1<0*)
g1=(4 (-256700232+1296384900 e2-970164375 (e2^2)-2198218750 (e2^3)+1543750000 (e2^4)+1406250000 (e2^5)))/(625 e2 (29488131-8127900 e2-113832500 (e2^2)-149000000 (e2^3)+400000000 (e2^4)));
g2=-((2 (-2134161091674039218421625601149503+140355978576685509466108656363508650 e2-1754283565687266330660162078035242500 (e2^2)+8737367532700018023712655697839625000 (e2^3)-21589750138203633706431942941437500000 (e2^4)+75104078381422305038823485272500000000 (e2^5)-513431090081409858991186790953125000000 (e2^6)+2056162232816920023362338942968750000000 (e2^7)-3838597863908140206223129570312500000000 (e2^8)+1617903538459282225849587890625000000000 (e2^9)+4615310325659104860692578125000000000000 (e2^10)-1399763372530238358330078125000000000000 (e2^11)-20086555324348486795898437500000000000000 (e2^12)+40652986288855743115234375000000000000000 (e2^13)-41911717827629003906250000000000000000000 (e2^14)+37465861770189697265625000000000000000000 (e2^15)-37513628819775390625000000000000000000000 (e2^16)+29057543024902343750000000000000000000000 (e2^17)-12113875488281250000000000000000000000000 (e2^18)+1593969726562500000000000000000000000000 (e2^19)+263671875000000000000000000000000000000 (e2^20)))/(1875 e2 (-3549+850 e2+10000 (e2^2)) (((25299+450 e2-295000 (e2^2)+250000 (e2^3)))^2) (-25299-450 e2-320000 (e2^2)+500000 (e2^3)) (29488131-8127900 e2-113832500 (e2^2)-149000000 (e2^3)+400000000 (e2^4)) (38082789-250500 e2-60467500 (e2^2)-515000000 (e2^3)+800000000 (e2^4))));
{c1,d1}
{c1,d1}//N
{g1,g2}//N

sol=Solve[{xd==0,yd==0},{x,y},Reals]//N

ClearAll[sol];
sol=Solve[{xd==0,yd==0},{x,y},Reals][[1]]//N
D[{xd,yd},{{x,y}}]/.sol;
Eigenvalues[D[{xd,yd},{{x,y}}]/.sol]//N

3. Perturbing c1 and e2 such that g1<0,g2>0, trace(J)>0

Quit

ClearAll[xd,yd,c1,d1,d2,e1,e2,f1,f2,g1,g2,sol];
xd=(x^2) y-x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y+x y+c1 (x^2)+d2 x-e2 y+f2;
f1=(1/2);f2=(1/10);
e1=(18/100);d2=(1/100);
d1=((c1 (e1-e2)+e2 f1+e1 (d2+f2))/e2);
c1=(233/1000);(*perturbed parameter for trace(J)>0*)
e2=(278/1000);(*perturbed parameter for g1<0*)
g1=(4 (-256700232+1296384900 e2-970164375 (e2^2)-2198218750 (e2^3)+1543750000 (e2^4)+1406250000 (e2^5)))/(625 e2 (29488131-8127900 e2-113832500 (e2^2)-149000000 (e2^3)+400000000 (e2^4)));
g2=-((2 (-2134161091674039218421625601149503+140355978576685509466108656363508650 e2-1754283565687266330660162078035242500 (e2^2)+8737367532700018023712655697839625000 (e2^3)-21589750138203633706431942941437500000 (e2^4)+75104078381422305038823485272500000000 (e2^5)-513431090081409858991186790953125000000 (e2^6)+2056162232816920023362338942968750000000 (e2^7)-3838597863908140206223129570312500000000 (e2^8)+1617903538459282225849587890625000000000 (e2^9)+4615310325659104860692578125000000000000 (e2^10)-1399763372530238358330078125000000000000 (e2^11)-20086555324348486795898437500000000000000 (e2^12)+40652986288855743115234375000000000000000 (e2^13)-41911717827629003906250000000000000000000 (e2^14)+37465861770189697265625000000000000000000 (e2^15)-37513628819775390625000000000000000000000 (e2^16)+29057543024902343750000000000000000000000 (e2^17)-12113875488281250000000000000000000000000 (e2^18)+1593969726562500000000000000000000000000 (e2^19)+263671875000000000000000000000000000000 (e2^20)))/(1875 e2 (-3549+850 e2+10000 (e2^2)) (((25299+450 e2-295000 (e2^2)+250000 (e2^3)))^2) (-25299-450 e2-320000 (e2^2)+500000 (e2^3)) (29488131-8127900 e2-113832500 (e2^2)-149000000 (e2^3)+400000000 (e2^4)) (38082789-250500 e2-60467500 (e2^2)-515000000 (e2^3)+800000000 (e2^4))));
{c1,d1}
{c1,d1}//N
{g1,g2}//N

Solve[{xd==0,yd==0},{x,y},Reals]//N

ClearAll[sol];
sol=Solve[{xd==0,yd==0},{x,y},Reals][[1]]//N
D[{xd,yd},{{x,y}}]/.sol;
Eigenvalues[D[{xd,yd},{{x,y}}]/.sol]//N

0.0003633093525180382` *2

Plotting the limit cycles when  g1<0,g2>0, trace(J)>0

Preparations

Quit

SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ParametricPlot,ListPlot,ListLinePlot};
SetDirectory[NotebookDirectory[]];
SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ListPlot,ListLinePlot,ListLogLogPlot,ParametricPlot,DateListPlot,DiscretePlot};
LaunchKernels[];

The function creating the plots

ClearAll[p,q,x,y,c1,d1,d2,e1,e2,f1,f2];
p[x_,y_]:=(x^2) y-x y-c1 (x^2)-d1 x+e1 y+f1;
q[x_,y_]:=-(x^2)y+x y+c1 (x^2)+d2 x-e2 y+f2;
f1=(1/2);f2=(1/10);
e1=(18/100);d2=(1/100);
d1=((c1 (e1-e2)+e2 f1+e1 (d2+f2))/e2);
c1=(233/1000);(*perturbed parameter for trace(J)>0*)
e2=(278/1000);(*perturbed parameter for g1<0*)
g1=(4 (-256700232+1296384900 e2-970164375 (e2^2)-2198218750 (e2^3)+1543750000 (e2^4)+1406250000 (e2^5)))/(625 e2 (29488131-8127900 e2-113832500 (e2^2)-149000000 (e2^3)+400000000 (e2^4)));
g2=-((2 (-2134161091674039218421625601149503+140355978576685509466108656363508650 e2-1754283565687266330660162078035242500 (e2^2)+8737367532700018023712655697839625000 (e2^3)-21589750138203633706431942941437500000 (e2^4)+75104078381422305038823485272500000000 (e2^5)-513431090081409858991186790953125000000 (e2^6)+2056162232816920023362338942968750000000 (e2^7)-3838597863908140206223129570312500000000 (e2^8)+1617903538459282225849587890625000000000 (e2^9)+4615310325659104860692578125000000000000 (e2^10)-1399763372530238358330078125000000000000 (e2^11)-20086555324348486795898437500000000000000 (e2^12)+40652986288855743115234375000000000000000 (e2^13)-41911717827629003906250000000000000000000 (e2^14)+37465861770189697265625000000000000000000 (e2^15)-37513628819775390625000000000000000000000 (e2^16)+29057543024902343750000000000000000000000 (e2^17)-12113875488281250000000000000000000000000 (e2^18)+1593969726562500000000000000000000000000 (e2^19)+263671875000000000000000000000000000000 (e2^20)))/(1875 e2 (-3549+850 e2+10000 (e2^2)) (((25299+450 e2-295000 (e2^2)+250000 (e2^3)))^2) (-25299-450 e2-320000 (e2^2)+500000 (e2^3)) (29488131-8127900 e2-113832500 (e2^2)-149000000 (e2^3)+400000000 (e2^4)) (38082789-250500 e2-60467500 (e2^2)-515000000 (e2^3)+800000000 (e2^4))));
{g1,g2}//N
{e2,e1,d2,d1,f1,f2,c1}
{e2,e1,d2,d1,f1,f2,c1}//N

ClearAll[nsol,nsol1,nsol2,ev,ev1,ev2,plotter];nsol=NSolve[Join@@Thread/@{{p[x,y],q[x,y]}==0,{x,y}>0},{x,y},20][[3]]
ev=Eigenvalues[D[{p[x,y],q[x,y]},{{x,y}}]/.nsol]
nsol1=NSolve[Join@@Thread/@{{p[x,y],q[x,y]}==0,{x,y}>0},{x,y},20][[1]]
ev1=Eigenvalues[D[{p[x,y],q[x,y]},{{x,y}}]/.nsol1]
nsol2=NSolve[Join@@Thread/@{{p[x,y],q[x,y]}==0,{x,y}>0},{x,y},20][[2]]
ev2=Eigenvalues[D[{p[x,y],q[x,y]},{{x,y}}]/.nsol2]
plotter[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:1000,ar_:Automatic,opts___]:=Module[{startingpoint,sys,solution,plot1},
startingpoint=({x,y}/.nsol)+shift;
sys:=NDSolveValue[Join[{u'[t]==p[u[t],v[t]],v'[t]==q[u[t],v[t]]},Thread[{u[0],v[0]}==startingpoint]],
{u,v},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
solution[t_]:=Through[sys[t]];
{ParametricPlot[Evaluate[solution[t]],{t,0,τ},Epilog->{Red,PointSize[0.05],Point[startingpoint],Orange,Point[{x,y}/.nsol],Blue,Point[{x,y}/.nsol1],Green,Point[{x,y}/.nsol2]},PlotRange->All,PlotPoints->pp,AspectRatio->ar,AxesLabel->{x,y},LabelStyle->Directive[14],ImageSize->200],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,x},LabelStyle->Directive[12],ImageSize->200],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,y},LabelStyle->Directive[12],ImageSize->200]}]

Plotter with arrow

ClearAll[nsol,nsol1,nsol2,ev,ev1,ev2,plotterarrow];nsol=NSolve[Join@@Thread/@{{p[x,y],q[x,y]}==0,{x,y}>0},{x,y},20][[3]]
ev=Eigenvalues[D[{p[x,y],q[x,y]},{{x,y}}]/.nsol]
nsol1=NSolve[Join@@Thread/@{{p[x,y],q[x,y]}==0,{x,y}>0},{x,y},20][[1]]
ev1=Eigenvalues[D[{p[x,y],q[x,y]},{{x,y}}]/.nsol1]
nsol2=NSolve[Join@@Thread/@{{p[x,y],q[x,y]}==0,{x,y}>0},{x,y},20][[2]]
ev2=Eigenvalues[D[{p[x,y],q[x,y]},{{x,y}}]/.nsol2]
plotterarrow[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:1000,ar_:Automatic,arrow_,opts___]:=Module[{startingpoint,sys,solution,plot1},
startingpoint=({x,y}/.nsol)+shift;
sys:=NDSolveValue[Join[{u'[t]==p[u[t],v[t]],v'[t]==q[u[t],v[t]]},Thread[{u[0],v[0]}==startingpoint]],
{u,v},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
solution[t_]:=Through[sys[t]];
{ParametricPlot[Evaluate[solution[t]],{t,0,τ},Epilog->{Black,Arrowheads->0.07,Arrow[{startingpoint,arrow}],Red,PointSize[0.05],Point[startingpoint],Orange,Point[{x,y}/.nsol],Blue,Point[{x,y}/.nsol1],Green,Point[{x,y}/.nsol2]},PlotRange->All,PlotPoints->pp,AspectRatio->ar,AxesLabel->{x,y},LabelStyle->Directive[14],ImageSize->200],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,x},LabelStyle->Directive[12],ImageSize->200],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,y},LabelStyle->Directive[12],ImageSize->200]}]

Figures

plotter[800,{0,0.47},Automatic,100,1000,1,Method->"BDF"]
plotterarrow[800,{0,0.47},Automatic,100,1000,1,{1.17,1.405},Method->"BDF"]

plotter[800,{0,0.46},Automatic,100,1000,1,Method->"BDF"]
plotterarrow[800,{0,0.46},Automatic,100,1000,1,{1.117,1.506},Method->"BDF"]

plotter[800,{0,0.05},Automatic,100,1000,1,Method->"BDF"]
plotterarrow[800,{0,0.05},Automatic,100,1000,1,{1.013,1.26},Method->"BDF"]

Model 2 with a stable outer an unstable inner limit cycle

The singular point is shifted into the origin

Singular points if x0=1

Quit

ClearAll[xd,yd,x,y,c1,d1,d2,e1,e2,f1,f2];
xd=(x^2) y-x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y+x y+c1 (x^2)+d2 x-e2 y+f2;
Solve[{xd==0,yd==0}/.x->1,{d1,y}]//FullSimplify

The singular point (if x0=1) is shifted into (0,0)

ClearAll[xd,yd,x,y,c1,d1,d2,e1,e2,f1,f2,x0,y0,x1d,y1d,xx1,yy1];
xd=(x^2) y-x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y+x y+c1 (x^2)+d2 x-e2 y+f2;
x0=1;
d1=((c1 (e1-e2)+e2 f1+e1 (d2+f2))/e2);y0=((c1+d2+f2)/e2);
xx1=x-x0;yy1=y-y0;
x1d=D[xx1,x]xd+D[xx1,y]yd/.{x->x1+x0,y->y1+y0}//Factor
y1d=D[yy1,x]xd+D[yy1,y]yd/.{x->x1+x0,y->y1+y0}//Factor

The Jacobian at the origin

Jac=D[{x1d,y1d},{{x1,y1}}];
JacOrigin=Jac/.{x1->0,y1->0}//Simplify

trace=Tr[JacOrigin]//Factor
Solve[trace==0,c1]//FullSimplify

c1=((d2-d2 e1-e2 (e2+f1)+f2-e1 f2)/(-1+e1+e2));
evalues=Eigenvalues[JacOrigin]//FullSimplify

beta=-((evalues[[1]])^2)//Factor

pp=x1d/.{x1->x,y1->y}//Factor
qq=y1d/.{x1->x,y1->y}//Factor

Lyapunov’s theorem

System

Quit

pp=(1/(-1+e1+e2))(-e2 x+e1 e2 x+(e2^2) x+d2 e1 (x^2)-e2 (x^2)+(e2^2) (x^2)-f1 (x^2)+e2 f1 (x^2)+e1 f2 (x^2)-e1 y+(e1^2) y+e1 e2 y-x y+e1 x y+e2 x y-(x^2) y+e1 (x^2) y+e2 (x^2) y);
qq=-(1/(-1+e1+e2))(d2 e1 x-e2 x-d2 e2 x+2 (e2^2) x-f1 x+2 e2 f1 x-f2 x+2 e1 f2 x+d2 e1 (x^2)-e2 (x^2)+(e2^2) (x^2)-f1 (x^2)+e2 f1 (x^2)+e1 f2 (x^2)-e2 y+e1 e2 y+(e2^2) y-x y+e1 x y+e2 x y-(x^2) y+e1 (x^2) y+e2 (x^2) y);

Program

Ser[s_]:=Plus@@Table[(x^i)(y^(s-i))p[i,s-i],{i,0,s}];
Hom[s_]:=Table[p[s-i,i],{i,0,s}];
hh=Sum[Ser[i],{i,2,6}];(*9*)
Lie=D[hh,x] pp+D[hh,y] qq//Expand;
RHS=g1((((x^2)+(y^2)))^2)+g2((((x^2)+(y^2)))^3)+g3((((x^2)+(y^2)))^4)//Expand;
vv=Lie-RHS//Expand;
CoefPol[f_,s_]:=Module[{m,lis,t},lis={}; m=Expand[f];  Do[Do[If[i+j==s,lis=AppendTo[lis,Coefficient[m,x^i y^j ]/.{x->0,y->0,z->0}]],{i,0,s}],{j,0,s}];
ls[s]=lis];
Do[CoefPol[vv,i],{i,1,9}];

Degree 1, 2

ls[1]
ls[2]//Factor;
sol2=Solve[ls[2]==0,Hom[2]]//Simplify;
{p[2,0],p[1,1],p[0,2]}={p[2,0],p[1,1],p[0,2]}/.sol2[[1]];
ls[2]//Simplify

Quadratic form

ClearAll[qv,mat,a11,det,eg];
qv=Ser[2]//FullSimplify
mat=1/2 D[qv,{{x,y},2}];
mat//MatrixForm
a11=mat[[1,1]]
det=Det[mat]//Factor
eg=Eigenvalues[mat]//Simplify//Factor;

Conditions for a positive definite quadratic form

d1=((c1 (e1-e2)+e2 f1+e1 (d2+f2))/e2);y0=((c1+d2+f2)/e2);
c1=((d2-d2 e1-e2 (e2+f1)+f2-e1 f2)/(-1+e1+e2));
beta=((d2 (e1^2)-e1 e2-d2 e1 e2+(e2^2)+e1 (e2^2)-(e2^3)-e1 f1+2 e1 e2 f1-e1 f2+2 (e1^2) f2)/(-1+e1+e2));

Reduce[ a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0&&f1>0 &&f2>0,{d2,e1,e2,f1,f2} ]//FullSimplify

Setting f1 and f2

ClearAll[f1,f2];
p[1,1]=1;f1=1;f2=1;

Reduce[ a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e1>0&&e2>0,{d2,e1,e2} ]//FullSimplify

Degree 3

ls[3]//Factor;
sol3=Solve[ls[3]==0,Hom[3]]//Simplify;
{p[3,0],p[2,1],p[1,2],p[0,3]}={p[3,0],p[2,1],p[1,2],p[0,3]}/.sol3[[1]]//Factor;
ls[3]//Simplify

Degree 4

ls[4]//Factor;
sol4=Solve[ls[4]==0,AppendTo[Hom[4],g1]];
{g1,p[4,0],p[3,1],p[2,2],p[1,3],p[0,4]}={g1,p[4,0],p[3,1],p[2,2],p[1,3],p[0,4]}/.sol4[[1]];
g1//Factor

Variables[g1]

Degree 5

ls[5]//Factor;
sol5=Solve[ls[5]==0,Hom[5]];
{p[5,0],p[4,1],p[3,2],p[2,3],p[1,4],p[0,5]}={p[5,0],p[4,1],p[3,2],p[2,3],p[1,4],p[0,5]}/.sol5[[1]]//Factor;
ls[5]//Factor

Degree 6

ls[6]//Factor;
sol6=Solve[ls[6]==0,AppendTo[Hom[6],g2]];
{p[6,0],p[5,1],p[4,2],p[3,3],p[2,4],p[1,5],p[0,6],g2}={p[6,0],p[5,1],p[4,2],p[3,3],p[2,4],p[1,5],p[0,6],g2}/.sol6[[1]]//Factor;

g2
Variables[g2]

p[2,2]=0;
g2;

Variables[g2]

Investigating the possible values of e1 and d2

ClearAll[d2,e1,e2];e1=(1/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2==0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify

ClearAll[d2,e1,e2];e1=(2/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2==0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify

ClearAll[d2,e1,e2];e1=(3/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify
Reduce[ g1==0&&g2==0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&d2>0&&e2>0,{d2,e2} ]//FullSimplify

Setting the value of d2

ClearAll[d2,e1,e2];d2=3/10;
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&e1≥0&&e2>0,{e1,e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&e1≥0&&e2>0,{e1,e2} ]//FullSimplify
Reduce[ g1==0&&g2==0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&e1≥0&&e2>0,{e1,e2} ]//FullSimplify

Setting the value of e1, d2 such that g1=0, g2<0

ClearAll[d2,e1,e2];e1=(3/10);d2=(4/10);
Reduce[ g1==0&&g2<0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&e2>0,{e2} ]//FullSimplify
Reduce[ g1==0&&g2>0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&e2>0,{e2} ]//FullSimplify
Reduce[ g1==0&&g2==0&&a11>0 && det>0&&beta>0&&y0>0&& c1>0&&d1>0&&e2>0,{e2} ]//FullSimplify

g1//Factor

g2//Factor

e2=-((NumericalApproximation[Root, (RowDefault["-18.9", FEExpressionsNumericalApproximationElider]), Root[13258566+1163754 #1+14544 ((#1)^2)-429 ((#1)^3)+24 ((#1)^4)+((#1)^5)&,1], ])/30);
e2//N

Perturbation

1. Setting e2 such that g1=0,g2<0, trace(J)=0

Quit

ClearAll[xd,yd,c1,d1,d2,e1,e2,f1,f2,g1,g2,sol];
xd=(x^2) y-x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y+x y+c1 (x^2)+d2 x-e2 y+f2;
f1=1;f2=1;
e1=(3/10);d2=(4/10);
d1=((c1 (e1-e2)+e2 f1+e1 (d2+f2))/e2);
c1=((d2-d2 e1-e2 (e2+f1)+f2-e1 f2)/(-1+e1+e2));
e2=-((NumericalApproximation[Root, (RowDefault["-18.9", FEExpressionsNumericalApproximationElider]), Root[13258566+1163754 #1+14544 ((#1)^2)-429 ((#1)^3)+24 ((#1)^4)+((#1)^5)&,1], ])/30);
g1=((-81843+215510 e2-80800 (e2^2)-71500 (e2^3)-120000 (e2^4)+150000 (e2^5))/(25 e2 (55851-60060 e2-125300 (e2^2)+28000 (e2^3)+160000 (e2^4))));
g2=-((-7163477245732849632+108834719548271657265 e2-618997960691032506300 (e2^2)+1654315932275893272000 (e2^3)-1921827260759476942500 (e2^4)+895027619630747687500 (e2^5)-5443787618252766375000 (e2^6)+23247943465166542500000 (e2^7)-40625221075278218750000 (e2^8)+34972974592196812500000 (e2^9)-20364241867282500000000 (e2^10)+26751631674075000000000 (e2^11)-36525826237125000000000 (e2^12)+20523696210000000000000 (e2^13)-5042947712500000000000 (e2^14)+20423760000000000000000 (e2^15)-40658469375000000000000 (e2^16)+34192681250000000000000 (e2^17)-13181250000000000000000 (e2^18)+1170000000000000000000 (e2^19)+450000000000000000000 (e2^20))/(300 e2 (-149+90 e2+200 (e2^2)) (((96-45 e2-325 (e2^2)+250 (e2^3)))^2) (-96+45 e2-200 (e2^2)+500 (e2^3)) (55851-60060 e2-125300 (e2^2)+28000 (e2^3)+160000 (e2^4)) (78749-72900 e2-169900 (e2^2)-60000 (e2^3)+320000 (e2^4))));
{c1,d1,e2}//N
{g1,g2}//N

Solve[{xd==0,yd==0},{x,y},Reals]//N

2. Perturbing e2 such that g1>0,g2<0, trace(J)=0

Quit

ClearAll[xd,yd,c1,d1,d2,e1,e2,f1,f2,g1,g2,sol];
xd=(x^2) y-x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y+x y+c1 (x^2)+d2 x-e2 y+f2;
f1=1;f2=1;
e1=(3/10);d2=(4/10);
d1=((c1 (e1-e2)+e2 f1+e1 (d2+f2))/e2);
c1=((d2-d2 e1-e2 (e2+f1)+f2-e1 f2)/(-1+e1+e2));
e2=(633/1000);(*perturbed parameter for g1>0*)
g1=((-81843+215510 e2-80800 (e2^2)-71500 (e2^3)-120000 (e2^4)+150000 (e2^5))/(25 e2 (55851-60060 e2-125300 (e2^2)+28000 (e2^3)+160000 (e2^4))));
g2=-((-7163477245732849632+108834719548271657265 e2-618997960691032506300 (e2^2)+1654315932275893272000 (e2^3)-1921827260759476942500 (e2^4)+895027619630747687500 (e2^5)-5443787618252766375000 (e2^6)+23247943465166542500000 (e2^7)-40625221075278218750000 (e2^8)+34972974592196812500000 (e2^9)-20364241867282500000000 (e2^10)+26751631674075000000000 (e2^11)-36525826237125000000000 (e2^12)+20523696210000000000000 (e2^13)-5042947712500000000000 (e2^14)+20423760000000000000000 (e2^15)-40658469375000000000000 (e2^16)+34192681250000000000000 (e2^17)-13181250000000000000000 (e2^18)+1170000000000000000000 (e2^19)+450000000000000000000 (e2^20))/(300 e2 (-149+90 e2+200 (e2^2)) (((96-45 e2-325 (e2^2)+250 (e2^3)))^2) (-96+45 e2-200 (e2^2)+500 (e2^3)) (55851-60060 e2-125300 (e2^2)+28000 (e2^3)+160000 (e2^4)) (78749-72900 e2-169900 (e2^2)-60000 (e2^3)+320000 (e2^4))));
{c1,d1}
{c1,d1}//N
{g1,g2}//N

sol=Solve[{xd==0,yd==0},{x,y},Reals]//N

3. Perturbing c1 and e2 such that g1>0,g2<0, trace(J)<0

Quit

ClearAll[xd,yd,c1,d1,d2,e1,e2,f1,f2,g1,g2,sol];
xd=(x^2) y-x y-c1 (x^2)-d1 x+e1 y+f1;
yd=-(x^2)y+x y+c1 (x^2)+d2 x-e2 y+f2;
f1=1;f2=1;
e1=(3/10);d2=(4/10);
d1=((c1 (e1-e2)+e2 f1+e1 (d2+f2))/e2);
c1=(79/100);(*perturbed parameter for trace(J)<0*)
e2=(633/1000);(*perturbed parameter for g1>0*)
g1=((-81843+215510 e2-80800 (e2^2)-71500 (e2^3)-120000 (e2^4)+150000 (e2^5))/(25 e2 (55851-60060 e2-125300 (e2^2)+28000 (e2^3)+160000 (e2^4))));
g2=-((-7163477245732849632+108834719548271657265 e2-618997960691032506300 (e2^2)+1654315932275893272000 (e2^3)-1921827260759476942500 (e2^4)+895027619630747687500 (e2^5)-5443787618252766375000 (e2^6)+23247943465166542500000 (e2^7)-40625221075278218750000 (e2^8)+34972974592196812500000 (e2^9)-20364241867282500000000 (e2^10)+26751631674075000000000 (e2^11)-36525826237125000000000 (e2^12)+20523696210000000000000 (e2^13)-5042947712500000000000 (e2^14)+20423760000000000000000 (e2^15)-40658469375000000000000 (e2^16)+34192681250000000000000 (e2^17)-13181250000000000000000 (e2^18)+1170000000000000000000 (e2^19)+450000000000000000000 (e2^20))/(300 e2 (-149+90 e2+200 (e2^2)) (((96-45 e2-325 (e2^2)+250 (e2^3)))^2) (-96+45 e2-200 (e2^2)+500 (e2^3)) (55851-60060 e2-125300 (e2^2)+28000 (e2^3)+160000 (e2^4)) (78749-72900 e2-169900 (e2^2)-60000 (e2^3)+320000 (e2^4))));
{c1,d1}
{c1,d1}//N
{g1,g2}//N

ClearAll[sol];
sol=Solve[{xd==0,yd==0},{x,y},Reals][[1]]//N
D[{xd,yd},{{x,y}}]/.sol;
Eigenvalues[D[{xd,yd},{{x,y}}]/.sol]//N

-0.0005995260663507829`*2

Plotting the limit cycles when  g1<0,g2>0, trace(J)>0

Preparations

Quit

SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ParametricPlot,ListPlot,ListLinePlot};
SetDirectory[NotebookDirectory[]];
SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ListPlot,ListLinePlot,ListLogLogPlot,ParametricPlot,DateListPlot,DiscretePlot};
LaunchKernels[];

The function creating the plots

ClearAll[p,q,x,y,c1,d1,d2,e1,e2,f1,f2];
p[x_,y_]:=(x^2) y-x y-c1 (x^2)-d1 x+e1 y+f1;
q[x_,y_]:=-(x^2)y+x y+c1 (x^2)+d2 x-e2 y+f2;
f1=1;f2=1;
e1=(3/10);d2=(4/10);
d1=((c1 (e1-e2)+e2 f1+e1 (d2+f2))/e2);
c1=(79/100);(*perturbed parameter for trace(J)<0*)
e2=(633/1000);(*perturbed parameter for g1>0*)
g1=((-81843+215510 e2-80800 (e2^2)-71500 (e2^3)-120000 (e2^4)+150000 (e2^5))/(25 e2 (55851-60060 e2-125300 (e2^2)+28000 (e2^3)+160000 (e2^4))));
g2=-((-7163477245732849632+108834719548271657265 e2-618997960691032506300 (e2^2)+1654315932275893272000 (e2^3)-1921827260759476942500 (e2^4)+895027619630747687500 (e2^5)-5443787618252766375000 (e2^6)+23247943465166542500000 (e2^7)-40625221075278218750000 (e2^8)+34972974592196812500000 (e2^9)-20364241867282500000000 (e2^10)+26751631674075000000000 (e2^11)-36525826237125000000000 (e2^12)+20523696210000000000000 (e2^13)-5042947712500000000000 (e2^14)+20423760000000000000000 (e2^15)-40658469375000000000000 (e2^16)+34192681250000000000000 (e2^17)-13181250000000000000000 (e2^18)+1170000000000000000000 (e2^19)+450000000000000000000 (e2^20))/(300 e2 (-149+90 e2+200 (e2^2)) (((96-45 e2-325 (e2^2)+250 (e2^3)))^2) (-96+45 e2-200 (e2^2)+500 (e2^3)) (55851-60060 e2-125300 (e2^2)+28000 (e2^3)+160000 (e2^4)) (78749-72900 e2-169900 (e2^2)-60000 (e2^3)+320000 (e2^4))));
{g1,g2}//N
{e2,e1,d2,d1,f1,f2,c1}
{e2,e1,d2,d1,f1,f2,c1}//N

ClearAll[nsol,ev,plotter];nsol=NSolve[Join@@Thread/@{{p[x,y],q[x,y]}==0,{x,y}>0},{x,y},20][[1]]
ev=Eigenvalues[D[{p[x,y],q[x,y]},{{x,y}}]/.nsol]
plotter[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:1000,ar_:Automatic,opts___]:=Module[{startingpoint,sys,solution,plot1},
startingpoint=({x,y}/.nsol)+shift;
sys:=NDSolveValue[Join[{u'[t]==p[u[t],v[t]],v'[t]==q[u[t],v[t]]},Thread[{u[0],v[0]}==startingpoint]],
{u,v},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
trafo[point_]:=point;
solution[t_]:=trafo[Through[sys[t]]];
{ParametricPlot[Evaluate[solution[t]],{t,0,τ},Epilog->{Red,PointSize[0.05],Point[startingpoint],Orange,Point[trafo[{x,y}/.nsol]]},PlotRange->All,PlotPoints->pp,AspectRatio->ar,AxesLabel->{x,y},LabelStyle->Directive[14],ImageSize->350],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,x},LabelStyle->Directive[12],ImageSize->250],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,y},LabelStyle->Directive[12],ImageSize->250]}]

Plotter with arrow

ClearAll[nsol,ev,plotterarrow];nsol=First@NSolve[Join@@Thread/@{{p[x,y],q[x,y]}==0,{x,y}>0},{x,y},20]
ev=Eigenvalues[D[{p[x,y],q[x,y]},{{x,y}}]/.nsol]
plotterarrow[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:1000,ar_:Automatic,arrow_,opts___]:=Module[{startingpoint,sys,solution,plot1},
startingpoint=({x,y}/.nsol)+shift;
sys:=NDSolveValue[Join[{u'[t]==p[u[t],v[t]],v'[t]==q[u[t],v[t]]},Thread[{u[0],v[0]}==startingpoint]],
{u,v},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
solution[t_]:=Through[sys[t]];
{ParametricPlot[Evaluate[solution[t]],{t,0,τ},Epilog->{Black,Arrowheads->0.07,Arrow[{startingpoint,arrow}],Red,PointSize[0.05],Point[startingpoint],Orange,
 Point[{x,y}/.nsol]
},PlotRange->All,PlotPoints->pp,AspectRatio->ar,AxesLabel->{x,y},LabelStyle->Directive[14],ImageSize->350],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,x},LabelStyle->Directive[12],ImageSize->250],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,y},LabelStyle->Directive[12],ImageSize->250]}]

Figures

Figure 12

plotter[300,{0,0.5},Automatic,100,1000,1,Method->"BDF"]
plotterarrow[300,{0,0.5},Automatic,100,1000,1,{1.214,3.559},Method->"BDF"]

Figure 13

plotter[800,{0,0.05},Automatic,100,1000,1,Method->"BDF"]
plotterarrow[800,{0,0.05},Automatic,100,1000,1,{1.05,3.396},Method->"BDF"]

Figure 14

plotter[600,{0,0.01},Automatic,100,1000,1,Method->"BDF"]
plotterarrow[600,{0,0.01},Automatic,100,1000,1,{1.005,3.458},Method->"BDF"]