So I am trying to use the MethodOfLines package to solve some equations which are somewhat similar to 2D diffusion equations. As a starting point, I have modified the heat equation example from the documentation to work in 2D, here is my code:

```
module Diffusion2D
using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets
# Parameters, variables, and derivatives
@parameters t x y
@variables u(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
Dy = Differential(y)
Dyy = Differential(y)^2
# 1D PDE and boundary conditions
eq = Dt(u(t, x, y)) ~ 0.1*(Dxx(u(t, x, y)) + Dyy(u(t, x, y)))
bcs = [u(0, x, y) ~ sin(2*pi*x) * sin(2*pi*y),
u(t, 0, y) ~ 0.0,
u(t, 1, y) ~ 0.0,
u(t, x, 0) ~ 0.0,
u(t, x, 1) ~ 0.0]
# Space and time domains
domains = [t ∈ Interval(0.0, 0.5),
x ∈ Interval(0.0, 1.0),
y ∈ Interval(0.0, 1.0)]
# PDE system
@named pdesys = PDESystem(eq, bcs, domains,[t, x, y],[u(t, x, y)])
# Method of lines discretization
dx = 0.01
dy = 0.01
order = 2
discretization = MOLFiniteDifference([x => dx, y => dy],t, approx_order=order)
# Convert the PDE problem into an ODE problem
prob = discretize(pdesys, discretization)
# Solve ODE problem
@time sol = solve(prob, Tsit5(), saveat=0.01)
end
```

This code seems to work correctly and it gives the right solution, and once it’s been compiled, the solve step runs in ~140 seconds typically, which is within reason for the dx/dy I am using. The main issue is that the first time I run the code, the solve step typically takes ~6000 seconds to compile and run (it was even longer than this before I put it all in a module). This strikes me as an unusually large compile time, have I made a mistake which is causing this issue? Or otherwise is there ways I can optimize the compile time? If it’s relevant, I am running this code inside a Jupyter notebook. Thanks!