Suppose we are evaluating the stability of a periodic solution of a dynamic system: $$ \dot{x}=f(x,t) $$ where $x\in \mathbb{R}^n$, and $f$ being a smooth vector field periodic in $t$ with $f(x,t) = f(x,t+T)$.
Suppose we find one periodic solution with the initial condition $x_0$.
I want to evaluate the stability of the obtained periodic solution, and it is determined by the eigenvalues of the corresponding monodromy matrix.
Then define a map $P:x(0)\mapsto x(T) $ that maps an initial point to one period later. Since $x_0$ is a period-1 point, we have $P(x_0)=x_0$, then the stability of the periodic point can be determined by the eigenvalues of the linearized map $\partial P/\partial x$ at $x_0$.
My question is, can I evaluate the linearized map $\partial P/\partial x$ numerically with the definition of derivatives in terms of limits?
Here is an example of the matrix element in the first row and first column:
$$ \partial P_1/\partial x_1 = \frac{P_1(x_0+[\epsilon,0,0,...,0])-P_1(x_0)}{\epsilon} $$
where $\epsilon$ is set to be small enough, and $P_1(x)$ can be numerically integrated.
Is this formula valid for all vector fields $f$, and for regular methods for numerical integration, such as the Runge-Kutta's methods?
Or does this formula make no sense?
My second question is, is the numerically obtained matrix $\partial P/\partial x$ the monodromy matrix that is transverse to the flow?
Thank you in advance!