For my research, I perform many differential corrections processes. Often times, for these processes I want the transition matrix, \Phi(t_1, t_0), which is obtained by propagating the linearization of the system alongside the nonlinear differential equations such that the following approximation can be made:

where \mathbf{\delta x}(t) is the variation with respect to the reference solution at time t. The linear differential equations are obviously

where A(t) is the jacobian of the nonlinear sytem at time t.

Now to the actual question, I see that there is infrastructure in DifferentialEquations.jl for finding sensitivity of the state to the parameters, but is there any ergonomic way to find this transition matrix given the equations of motion and Jacobian without gluing them together in some function? My goal is being able to do this in general for many different systems which only need to provide their EOM and Jacobian definitions. A contrived small example would be

```
abstract type AbstractSystem end
struct SystemA <: AbstractSystem
alpha::Float64
end
struct SystemB <: AbstractSystem
beta::Float64
end
function eom!(du, u, p::SystemA, t)
du[1] = u[1] + exp(u[2])
du[2] = -p.alpha* u[2]
end
function jac!(J, u, p::SystemA, t)
J[1, 1] = 1
J[1, 2] = exp(u[2])
J[2, 1] = 0
J[2, 2] = -p.alpha
end
function eom!(du, u, p::SystemB, t)
du[1] = u[2] - u[2]^3 - p.beta* u[1]
du[2] = u[1] - u[2] - u[1] * u[2]
end
function jac!(J, u, p::SystemB, t)
J[1,1] = -p.beta
J[1,2] = 1 - 3 * u[2]^2
J[2,1] = 1 - u[2]
J[2,2] = -1 - u[1]
end
function stm(u0, system::T, tspan) where {T <: AbstractSystem}
# Do something here to get the `sol` as well as the
# transition matrices over timespan.
end
```