Return to Answer

For more details see: Andronov-Hopf bifurcation.

In case anyone wishes to collaborate on these topics, please, communicate.

For more details, see: Andronov-Hopf bifurcation.

MatrixForm@FullSimplify[F0[X, μ] - (JX0μ0.X + 1/2! BB[X, X] + 1/3! CC[X, X, X]) /. {x -> t x, y -> t y}]
(*{0,0}*)

(*Eigenvectors for J[X0,μ0]*)
vp = ComplexExpand@Eigenvectors[JX0μ0]
q = vp[[2]];
qc = vp[[1]];
MatrixForm@q
MatrixForm@qc
MatrixForm@Simplify[JX0μ0.q - I ω q]
(*Eigenvectors for Transpose[J[X0,μ0]]*)
vpt = ComplexExpand[Eigenvectors[JtX0μ0]]
(*Normalization constant*)
cn = ComplexExpand[Conjugate[vp[[1]] . vpt[[1]]]]
p = Expand@Simplify[vpt[[2]]/cn];
pc = ComplexExpand[Conjugate[p]];
MatrixForm@p
MatrixForm@pc
Simplify[JtX0μ0.p-I ω p]

Factor@ComplexExpand[Re[1/(2 ω) (p.CC[q, q, qc] - 2 (p.BB[q, Inverse[JX0μ0].BB[q, qc]]) + p.BB[qc, Inverse[2 I ω*IdentityMatrix[2] - JX0μ0].BB[q, q]])] /. ω -> ω0]
(*-((2 + α^2)/(2 α (1 + α^2)))*)

A=JX0μ0; (*linear approximation*)
MatrixForm@FullSimplify[F0[X, μ] - (A.X + 1/2! BB[X, X] + 1/3! CC[X, X, X]) /. {x -> t x, y -> t y}]
(*{0,0}*)

(*Eigenvectors of A=J[X0,μ0]*)
vp = ComplexExpand@Eigenvectors[A]
q = vp[[2]];
qc = vp[[1]];
MatrixForm@q
MatrixForm@qc
MatrixForm@Simplify[A.q - I ω q]
(*Eigenvectors of Transpose[A]*)
At=JtX0μ0;
vpt = ComplexExpand[Eigenvectors[At]]
(*Normalization constant*)
cn = ComplexExpand[Conjugate[vp[[1]] . vpt[[1]]]]
p = Expand@Simplify[vpt[[2]]/cn];
pc = ComplexExpand[Conjugate[p]];
MatrixForm@p
MatrixForm@pc
Simplify[At.p-I ω p]

Factor@ComplexExpand[Re[1/(2 ω) (p.CC[q, q, qc] - 2 (p.BB[q, Inverse[A].BB[q, qc]]) + p.BB[qc, Inverse[2 I ω*IdentityMatrix[2] - A].BB[q, q]])] /. ω -> ω0]
(*-((2 + α^2)/(2 α (1 + α^2)))*)

f1[x_, y_] := α - (β + 1) x + x^2 y;
f2[x_, y_] := β x - x^2 y;
F[{x_, y_}, {α_, β_}] := Evaluate@{f1[x, y], f2[x, y]};
X = {x, y};
μ = {α, β};

Factor@ComplexExpand[Re[1/(2 ω) (p.D3X0μ0.q.q.qc- 2(p.D2X0μ0.q.Inverse[JX0μ0].D2X0μ0.q. qc) +
 p.D2X0μ0.qc. Inverse[2 I ω*IdentityMatrix[2]-JX0μ0].D2X0μ0.q. q)] /. ω -> ω0]
(*-((2 + α^2)/(2 α (1 + α^2)))*)

f1[x_, y_] := α - (β + 1) x + x^2 y;
f2[x_, y_] := β x - x^2 y;
F[{x_, y_}, {α_, β_}] := Evaluate@{f1[x, y], f2[x, y]};
X = {x, y};
μ = {α, β};
U = {u, v};
R = {r, s};

Factor@ComplexExpand[Re[1/(2 ω) (p.CC[q, q, qc] - 2 (p.BB[q, Inverse[JX0μ0].BB[q, qc]]) + p.BB[qc, Inverse[2 I ω*IdentityMatrix[2] - JX0μ0].BB[q, q]])] /. ω -> ω0]
(*-((2 + α^2)/(2 α (1 + α^2)))*)