Consider an exemplary PDE with temporal & spatial derivatives

\partial_t u(t, x) + D_x f\big(u(t, x)\big) = 0

Where D_x is a general differential operator, like the Laplacian and f is some function.

I can use MethodOfLines.jl to obtain a semi-discretization U(t) = F\big(U(t)\big)
through invoking discretize(sys::PDESystem, disc::D) which gives a AbstractSciMLProblem. Is it possible to call a routine from one of the other SiML packages to obtain the Jacobian J(U) = \frac{\partial F\big(U(t)\big)}{\partial U} ?
Ideally, this would be possible via algorithmic/automated differentiation, although Finite DIfferences would also be a suitable starting point.

During discretize if you pass jac=true then it will generate it, and sys.jac[] or prob.f.jac are these pieces. @xtalax this seems to be missing in the tutorials?

Now I would like to evaluate that thing to obtain an actual matrix. I guess one of the arguments arg1, arg2 is the point of evaluation U^* and t is probably time t if the Jacobian is time-dependent, i.e., J\big(t, U(t) \big).
Question is what arg1, arg2 are - I will see if I can figure this out via the ModelingToolkit documentation.