Floquet theory for periodic delay differential equations: current numerical routines

Question

I would like to determine the stability of a system of periodic delay differential equations (a seasonal host-parasite model). I've tried to implement the method described in Lemma 2.5 in this paper:

• Liang, X., Zhang, L. & Zhao, XQ. Basic Reproduction Ratios for Periodic Abstract Functional Differential Equations (with Application to a Spatial Model for Lyme Disease). J Dyn Diff Equat 31, 1247–1278 (2019).:

Lemma 2.5 Assume that (E, E₊) is an ordered Banach space with E₊ being normal and Int(E₊) $\neq$ ∅, which is equipped with the norm ||·||_E . Let L be a positive bounded linear operator. Choose v₀ ∈ Int(E₊) and define a_n = ||Lv_n-1||_E , v_n = $\frac{Lv_{n-1}}{ > a_{n}}$, ∀n ≥ 1. If $\lim_{(n→+∞)}a_n$ exists, then r(L) = $\lim_{n→+∞}a_n$.

But I'm not sure it's working properly (I've probably made a mistake somewhere in the implementation). The problem is it doesn't always seem to correspond to equilibrium points I find when simply simulating the model for a very long time. I looked into different ways of approximating Floquet multipliers/monodromy matrices etc. for delay equations, just so I could double check my work, and there seem to be many different approaches (e.g., semi-discretization, finite element analysis, Chebyshev polynomials, etc.)

My question, for anyone much more familiar with this literature than I am, is: Is there a currently preferred method? Are there any that are easy and quick to implement (I'm not a very sophisticated programmer or mathematician)? I would be especially interested in an example with published code. Some advice would be greatly appreciated.