diff --git a/Project.toml b/Project.toml
index b7ba701..bebeb8c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.1.4"
 
 [deps]
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
+FunctionMaps = "a85aefff-f8ca-4649-a888-c8e5398bc76c"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
@@ -14,12 +15,13 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 
 [compat]
-DomainSets = "0.6, 0.7"
+DomainSets = "0.7.17, 0.8"
+FunctionMaps = "0.1"
 Interpolations = "0.14, 0.15, 0.16"
 Lux = "1"
 Markdown = "1"
-MethodOfLines = "0.11.12, 0.12"
-ModelingToolkit = "9.7.1, 10, 11.27"
+MethodOfLines = "0.11.14, 0.12"
+ModelingToolkit = "11.30.1"
 NonlinearSolve = "4.12"
 OptimizationOptimJL = "0.3, 0.4"
 OrdinaryDiffEq = "6.80.0, 7"
diff --git a/lib/burgers.jl b/lib/burgers.jl
index f7119c5..f20b1f6 100644
--- a/lib/burgers.jl
+++ b/lib/burgers.jl
@@ -70,71 +70,40 @@ function burgers_2d()
     u_exact(x, y, t) = 3 / 4 - 1 / (4 * (1 + exp(R * (-t - 4x + 4y) / 32)))
     v_exact(x, y, t) = 3 / 4 + 1 / (4 * (1 + exp(R * (-t - 4x + 4y) / 32)))
     analytic = [
-        u(x, y, t) ~ u_exact(x, y, t),
-        v(x, y, t) ~ v_exact(x, y, t),
+        u(t, x, y) ~ u_exact(x, y, t),
+        v(t, x, y) ~ v_exact(x, y, t),
     ]
 
     eq = [
-        Dt(u(x, y, t)) + u(x, y, t) * Dx(u(x, y, t)) + v(x, y, t) * Dy(u(x, y, t)) ~
-            (1 / R) *
-            (
-            Dxx(
-                u(
-                    x,
-                    y,
-                    t
-                )
-            ) +
-                Dyy(
-                u(
-                    x,
-                    y,
-                    t
-                )
-            )
-        ),
-        Dt(v(x, y, t)) + u(x, y, t) * Dx(v(x, y, t)) + v(x, y, t) * Dy(v(x, y, t)) ~
-            (1 / R) *
-            (
-            Dxx(
-                v(
-                    x,
-                    y,
-                    t
-                )
-            ) +
-                Dyy(
-                v(
-                    x,
-                    y,
-                    t
-                )
-            )
-        ),
+        Dt(u(t, x, y)) + u(t, x, y) * Dx(u(t, x, y)) + v(t, x, y) * Dy(u(t, x, y)) ~
+            (1 / R) * (Dxx(u(t, x, y)) + Dyy(u(t, x, y))),
+        Dt(v(t, x, y)) + u(t, x, y) * Dx(v(t, x, y)) + v(t, x, y) * Dy(v(t, x, y)) ~
+            (1 / R) * (Dxx(v(t, x, y)) + Dyy(v(t, x, y))),
     ]
 
     domains = [
+        t ∈ Interval(t_min, t_max),
         x ∈ Interval(x_min, x_max),
         y ∈ Interval(y_min, y_max),
-        t ∈ Interval(t_min, t_max),
     ]
 
     bcs = [
-        u(x, y, 0) ~ u_exact(x, y, 0),
-        u(0, y, t) ~ u_exact(0, y, t),
-        u(x, 0, t) ~ u_exact(x, 0, t),
-        u(1, y, t) ~ u_exact(1, y, t),
-        u(x, 1, t) ~ u_exact(x, 1, t), v(x, y, 0) ~ v_exact(x, y, 0),
-        v(0, y, t) ~ v_exact(0, y, t),
-        v(x, 0, t) ~ v_exact(x, 0, t),
-        v(1, y, t) ~ v_exact(1, y, t),
-        v(x, 1, t) ~ v_exact(x, 1, t),
+        u(0, x, y) ~ u_exact(x, y, 0),
+        u(t, 0, y) ~ u_exact(0, y, t),
+        u(t, x, 0) ~ u_exact(x, 0, t),
+        u(t, 1, y) ~ u_exact(1, y, t),
+        u(t, x, 1) ~ u_exact(x, 1, t),
+        v(0, x, y) ~ v_exact(x, y, 0),
+        v(t, 0, y) ~ v_exact(0, y, t),
+        v(t, x, 0) ~ v_exact(x, 0, t),
+        v(t, 1, y) ~ v_exact(1, y, t),
+        v(t, x, 1) ~ v_exact(x, 1, t),
     ]
 
     tags = ["2D", "Non-Monotonic", "Viscous", "Burgers", "Advection", "Dirichlet"]
 
     @named burgers_2d = PDESystem(
-        eq, bcs, domains, [t, x, y], [u(x, y, t), v(x, y, t)],
+        eq, bcs, domains, [t, x, y], [u(t, x, y), v(t, x, y)],
         analytic = analytic, metadata = tags
     )
 
diff --git a/lib/linear_diffusion.jl b/lib/linear_diffusion.jl
index b97b6ee..630c48f 100644
--- a/lib/linear_diffusion.jl
+++ b/lib/linear_diffusion.jl
@@ -12,10 +12,10 @@ The equation is given by:
 ```
 """
 function heat_1d1()
-    @variables x t u(..)
-    @parameters D
+    @parameters x t D
+    @variables u(..)
 
-    Dxx = Differential(x)
+    Dxx = Differential(x)^2
     Dt = Differential(t)
 
     eqs = [Dt(u(t, x)) ~ D * Dxx(u(t, x))]
diff --git a/src/PDESystemLibrary.jl b/src/PDESystemLibrary.jl
index fc413a9..2d1d3fa 100644
--- a/src/PDESystemLibrary.jl
+++ b/src/PDESystemLibrary.jl
@@ -1,5 +1,5 @@
 module PDESystemLibrary
-using ModelingToolkit, DomainSets
+using ModelingToolkit, DomainSets, FunctionMaps
 using OrdinaryDiffEq
 using OrdinaryDiffEqSDIRK
 using Interpolations
diff --git a/test/mol_test.jl b/test/mol_test.jl
index f4aff06..8655014 100644
--- a/test/mol_test.jl
+++ b/test/mol_test.jl
@@ -14,16 +14,10 @@ N = 100
 # them requires either an upwind/finite-volume discretization or an explicit
 # solver such as `Vern9()`.
 #
-# `:advdiff3` is also unstable, but only on Julia >= 1.13: on earlier versions
-# the integrator's adaptive step keeps it under control, while changes to
-# floating-point tie-breaking / lowering on 1.13 push the solver over the
-# stiffness threshold and FBDF aborts with `dt < eps`. Tracked conditionally
-# so we keep regression coverage on lts/1.x and don't false-fail on 1.13+.
-const BROKEN_EXAMPLES = if VERSION >= v"1.13-"
-    Set([:adv3, :advdiff3])
-else
-    Set([:adv3])
-end
+const BROKEN_EXAMPLES = Set([:adv3])
+# `:advdiff3` has diffusion, but the third-order term still makes the default
+# FBDF solve hit `dt < eps`; Rodas5P keeps this example covered.
+const EXAMPLE_ALGORITHMS = Dict(:advdiff3 => Rodas5P())
 
 for ex in PSL.all_systems
     @testset "Example: $(ex.name)" begin
@@ -63,7 +57,7 @@ for ex in PSL.all_systems
                     false
                 end
             else
-                sol = solve(prob, FBDF())
+                sol = solve(prob, get(EXAMPLE_ALGORITHMS, ex.name, FBDF()))
                 @test sol.retcode == SciMLBase.ReturnCode.Success
             end
         end