I am trying to numerically calculate the Lyapunov exponents for the seir epidemic model given as: $$ s^{'} = b - bs - \beta si \\ e^{'} = \beta si - (\alpha + b)e \\ i^{'} = \alpha e - (\gamma + b) i \\ r^{'} = \gamma i - br \\ s + e + i + r = 0 $$
here, $\beta$ is the control parameter with $\beta = \beta_0(1+\beta_1 cos2\pi t)$.
Initially, I compute the derivatives as $$ \frac{ds(t)}{ds} = - b - \beta i \\ \frac{de(t)}{de} = \beta s - (\alpha + b) \\ \frac{di(t)}{di} = \alpha - (\gamma + b) \\ \frac{dr(t)}{dr} = - b \\ $$
After which, to compute the Lyapunov exponents, I use the formula $\lambda = lim_{t \rightarrow \infty}\frac{1}{t} ln |\frac{dF(t)}{dF}|$, where $\frac{dF(t)}{dF}$ are taken as the four equations represented above. As a result, I get four Lyapunov exponents. I integrate the system using an RK4 integrator for a given time period and average to find the maximal Lyapunov exponents for a chosen parameter $\beta$. However, I am certain the values that I attain are wrong since I get two positive Lyapunov exponents and two negative Lyapunov exponents, even for a periodic attractor. I am also confused about calculating the derivative since e and i are dependent on each other.
I think I do not grasp the theory very well. Can someone help me understand this better?