Find Jacobian within ode-solver

Suppose I have an autonomous system of differential equations:

Equation 1: dx(t)/dt=f(x(t))

When I solve it in Matlab via one of the ode-solvers I get a discrete solution vector corresponding to discrete time steps. So far, so good.

Now I want to solve another system of differential equations:

Equation 2: dA(x(t),t)/dt=J(x(t),t)*A(x(t),t)

A and J are matrices. Now here is the catch: J(x(t),t) is the Jacobian of f(x(t)) evaluated at a certain x(t) and t. For example for the first step t is 0, so I need to find the Jacobian of f(x(t)) evaluated at x(0),0: J(x(0),0)

First off: I have tried to make a function including the actual (analytical) Jacobian. However my system is quite large, so the calculations are huge. Now I want to try finding the Jacobian numerically. I suppose I could try simple finite differences.

I was made aware of the numjac-function which - to my knowledge - computes the Jacobian of the right-hand-side of my system numerically.

Now to my problem: I have trouble at implementing. It seems to me that I can't solve equations 1 and 2 seperately, because I need the value of x(t) for a certain (t) [to find J(x(t),t)] and the only time I have this value is when Equation 2 uses the same time step as Equation 1. Thus I have tried merging the two equations into one (simply by adding to the vector x) and then finding the Jacobian J(x(t),t) via numjac inside the equation which is solved via an ode-solver. This works...somewhat. But it is still far too slow. I tried exporting it to C-Code (building a .mex-file of my equation via the coder) but this won't work because of the numjac-function inside.

I have the feeling my way is a bit circuitous...Is there a more elegant way to my problem? Maybe a way to solve both equations independently?

Sorry for the long post.