Are there public Mathematica scripts for computing Lyapunov's first focal value for planar systems?

Question

This seemed to me a utility which must have been created already hundreds of times, but I only found an implementation of 2022, which did not work on my example, so I needed to fix it first. Here's my example

ClearAll["Global\[OpenCurlyQuote]*"];
fg = {Subscript[f, 1] - Subscript[f, 1] x - (
   x (Subscript[d, 
       2] (2 + Subscript[e, 1]) + (2 + Subscript[e, 1]) Subscript[f, 
       2] + Subscript[c, 
       1] (Subscript[e, 1] - Subscript[e, 
         2] + (2 + Subscript[e, 2]) x)))/(
   2 + Subscript[e, 2]) + (Subscript[e, 1] + x + x^2) y, 
  Subscript[f, 2] + 
   x (Subscript[d, 2] + Subscript[c, 1] x) - (Subscript[e, 2] + x + 
      x^2) y}; var = {x, y};
equ = Solve[Thread[fg == 0], var][[1]] // FullSimplify
Print["Hopf necessary cond"]
J = D[fg, {var}] /. equ;
tr = Tr[J];
cc = Solve[tr == 0, Subscript[c, 1]] // Flatten // FullSimplify
fgc = fg /. cc;
Print["equilibrium"]
equ = Solve[Thread[fgc == 0], var][[1]] // FullSimplify
J = D[fgc, {var}] /. equ // FullSimplify;
det = Det[J];

My script below works now, after some help I got from the question 114769 Package functions and symbolic calculations. However, I will leave the question, with the hope that a better script than mine will be suggested by someone for this important object .

L1Planar[fg_, var_, equilibrium_ : {}] := 
  Module[{J, xyshift, Tm, Tinvuv, FG,
    derivatives, a, b, u, v, i, j, L1}, {x, y} = var;
   J = Simplify[D[fg, {var}] /. equilibrium];
   xyshift = {x -> x + (x /. equilibrium), 
     y -> y + (y /. equilibrium)};
   Tm = {{1, 0}, {-a/ome, -b/ome}};
   Tinvuv = Inverse[Tm] . {u, v};
   FG = (Tm . fg /. xyshift) 
      /. {x -> Tinvuv[[1]], y -> Tinvuv[[2]]} /. {a -> J[[1, 1]], 
      b -> J[[1, 2]]};
   derivatives = {};
   For[i = 0, i <= 3, i++, For[j = 0, j <= 3 - i, j++, 
     derivatives = 
      Join[derivatives, 
       {Subscript[F, i, 
          j] -> (D[FG[[1]], {u, i}, {v, j}] /. {u -> 0, v -> 0}), 
        Subscript[G, i, 
          j] -> (D[FG[[2]], {u, i}, {v, j}] /. {u -> 0, v -> 0})}]]];
   L1 = Subscript[F, 3, 0] + Subscript[F, 1, 2] + 
     Subscript[G, 0, 3] +
     Subscript[G, 2, 1] + 
     1/ome (Subscript[F, 1, 
          1] (Subscript[F, 2, 0] + Subscript[F, 0, 2]) - 
        Subscript[G, 1, 1] (Subscript[G, 2, 0] + Subscript[G, 0, 2]) +
         Subscript[F, 0, 2] Subscript[G, 0, 2] - 
        Subscript[F, 2, 0] Subscript[G, 2, 0]);
   L1 /. derivatives]; 
L1 = L1Planar[fgc, var, equ] // FullSimplify
L1d = L1 //. {ome^2 -> det, ome^-2 -> det^(-1)}