I am interested in analyzing the stability of the periodic orbits resulting from the Van der Pol system periodically perturbed by a time-dependent external forcing. Mathematically it would be the following system of nonlinear differential equations
$\begin{equation} \left \{ \begin{aligned} \frac{dx}{dt} & = \mu \left ( x - \frac{x^3}{3} - y \right ) + \varepsilon\cos(\omega t), \\ \frac{dy}{dt} & = \frac{x}{\mu}. \end{aligned} \right. \end{equation}$
When the system is autonomous, i.e $\varepsilon = 0$, I only have to solve the matrix differential equation $\dot{\Phi}(t) = D_{x}\mathbf{f}(X(t))\Phi(t)$ with the initial condition $\Phi(0) = \mathbf{I}$, where $D_{x}\mathbf{f}$ is the Jacobian matrix of the vector field that defines the system of differential equations. $\Phi(T)$ is the monodromy matrix of the periodic orbit $X(t)$, where $T$ represents its period. How can I calculate the monodromy matrix of the non-autonomous system?