MethodOfLines: How can I compute the Jacobian of the system of ODEs?

Specific Domains Modelling & Simulations
,
1

How can I compute the Jacobian of the system of ODEs obtained by spatially discretising a linear advection equation? A MWE follows:

using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets

@independent_variables t x
@variables T(..)
Dₜ = Differential(t)
Dₓ = Differential(x)

eq  = Dₜ(T(t, x)) ~ -Dₓ(T(t, x))

T₀ = 10.0
bcs = [T(0, x) ~ T₀,
       T(t, 0) ~ T₀+exp(-((t-1.0)^2)/0.1)]

domains = [t ∈ Interval(0.0, 1.0),
           x ∈ Interval(0.0, 2.0)]

@named pdesys = PDESystem(eq, bcs, domains, [t, x], [T(t, x)])

Δx = 0.1
discretization = MOLFiniteDifference([x => Δx], t, advection_scheme = UpwindScheme())

prob = discretize(pdesys,discretization)

How do I get the Jacobian matrix from the prob variable, which is of type ODEProblem?

An identical question has been asked and answered some two years ago: MethodOfLines.jl: Obtain Jacobian of MOL semi-discretization. The advise given by @ChrisRackauckas there was “During discretize if you pass jac=true then it will generate it, and sys.jac or prob.f.jac are these pieces.” But when I try

prob = discretize(pdesys,discretization,jac=true)

Julia seems to never stop computing. I confess that after some fifteen minutes I had to kill it. The more so that I cannot find such option anywhere in the docs. In fact, I cannot find any docs for the discretize function at all.

Combined with the fact that the advice is now some two years old (maybe it is now deprecated?), that is also why I decided not to “necropost” but rather write a new post.