diff --git a/Project.toml b/Project.toml
index 00ebfcb08..5b339a461 100644
--- a/Project.toml
+++ b/Project.toml
@@ -30,7 +30,7 @@ DomainSets = "0.7"
 IfElse = "0.1"
 Interpolations = "0.14, 0.15, 0.16"
 Latexify = "0.15, 0.16"
-ModelingToolkit = "9.17"
+ModelingToolkit = "10"
 OrdinaryDiffEq = "6"
 PDEBase = "0.1.17"
 PrecompileTools = "1"
diff --git a/docs/Project.toml b/docs/Project.toml
index b5062c4dc..55cecc549 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -33,7 +33,7 @@ IfElse = "0.1"
 Interpolations = "0.14, 0.15, 0.16"
 Latexify = "0.15, 0.16"
 MethodOfLines = "0.11"
-ModelingToolkit = "9"
+ModelingToolkit = "10"
 OrdinaryDiffEq = "6"
 PDEBase = "0.1.11"
 Plots = "1"
diff --git a/docs/src/howitworks.md b/docs/src/howitworks.md
index f1eb264be..2e4d1fb19 100644
--- a/docs/src/howitworks.md
+++ b/docs/src/howitworks.md
@@ -12,4 +12,4 @@ Next, the boundary conditions are discretized, creating an equation for each poi
 
 After that, the system of PDEs is discretized, creating a finite difference equation for each point in their interior. Specific terms are recognized, and the best implemented scheme for these terms dispatched. For example, advection terms are discretized with the upwind scheme. There are also special schemes for the nonlinear Laplacian and spherical Laplacian. See [here for how this term matching occurs](https://github.com/SciML/MethodOfLines.jl/blob/master/src/discretization/generate_finite_difference_rules.jl), note that the order the generated rules are applied is important, with more specific rules applied first to avoid their terms being matched incorrectly by more general rules.  The [SymbolicUtils.jl docs](https://symbolicutils.juliasymbolics.org/rewrite/) are a useful companion here. See [here for the practical implementation of the finite difference schemes](https://github.com/SciML/MethodOfLines.jl/blob/master/src/discretization/differential_discretizer.jl).
 
-Now we have a system of equations which are either ODEs, linear, or nonlinear equations and an equal number of unknowns. See [here for the system](@ref brusssys) that is generated for the Brusselator at low point count. The structure of the system is simplified with `ModelingToolkit.structural_simplify`, and then either an `ODEProblem` or `NonlinearProblem` is returned. Under the hood, the `ODEProblem` generates a fast semidiscretization, written in Julia with `RuntimeGeneratedFunctions`. See [here for an example of the generated code](@ref brusscode) for the Brusselator system at low point count.
+Now we have a system of equations which are either ODEs, linear, or nonlinear equations and an equal number of unknowns. See [here for the system](@ref brusssys) that is generated for the Brusselator at low point count. The structure of the system is simplified with `ModelingToolkit.mtkcompile`, and then either an `ODEProblem` or `NonlinearProblem` is returned. Under the hood, the `ODEProblem` generates a fast semidiscretization, written in Julia with `RuntimeGeneratedFunctions`. See [here for an example of the generated code](@ref brusscode) for the Brusselator system at low point count.
diff --git a/src/MOL_discretization.jl b/src/MOL_discretization.jl
index 79f179f4d..801b54a07 100644
--- a/src/MOL_discretization.jl
+++ b/src/MOL_discretization.jl
@@ -66,15 +66,15 @@ function get_discrete(pdesys, discretization)
         [Num(x) => s.grid[x] for x in s.x̄], [Num(u) => s.discvars[u] for u in s.ū]))
 end
 
-function ModelingToolkit.ODEFunctionExpr(
+function ODEFunctionExpr(
         pdesys::PDESystem, discretization::MethodOfLines.MOLFiniteDifference)
     sys, tspan = SciMLBase.symbolic_discretize(pdesys, discretization)
     try
         if tspan === nothing
             @assert true "Codegen for NonlinearSystems is not yet implemented."
         else
-            simpsys = structural_simplify(sys)
-            return ODEFunctionExpr(simpsys)
+            simpsys = mtkcompile(sys)
+            return ODEFunction(simpsys; expression = Val{true})
         end
     catch e
         println("The system of equations is:")
@@ -94,7 +94,7 @@ function SciMLBase.ODEFunction(
         if tspan === nothing
             @assert true "Codegen for NonlinearSystems is not yet implemented."
         else
-            simpsys = structural_simplify(sys)
+            simpsys = mtkcompile(sys)
             if analytic !== nothing
                 analytic = analytic isa Dict ? analytic : Dict(analytic)
                 s = getfield(sys, :metadata).discretespace
diff --git a/src/discretization/staggered_discretize.jl b/src/discretization/staggered_discretize.jl
index 1d4b80c19..8ff0c68fb 100644
--- a/src/discretization/staggered_discretize.jl
+++ b/src/discretization/staggered_discretize.jl
@@ -3,7 +3,7 @@ function SciMLBase.discretize(pdesys::PDESystem,
         analytic = nothing, kwargs...) where {G <: StaggeredGrid}
     sys, tspan = SciMLBase.symbolic_discretize(pdesys, discretization)
     try
-        simpsys = structural_simplify(sys)
+        simpsys = mtkcompile(sys)
         if tspan === nothing
             add_metadata!(get_metadata(sys), sys)
             return prob = NonlinearProblem(simpsys, ones(length(simpsys.unknowns));
diff --git a/test/pde_systems/MOL_NonlinearProblem.jl b/test/pde_systems/MOL_NonlinearProblem.jl
index bc219ad15..83c2b5c1c 100644
--- a/test/pde_systems/MOL_NonlinearProblem.jl
+++ b/test/pde_systems/MOL_NonlinearProblem.jl
@@ -14,7 +14,7 @@ using DomainSets
         0 ~ b + d - 2 * c,
         0 ~ d - 1]
     @named ns = NonlinearSystem(eqs, [a, b, c, d], [])
-    sns = structural_simplify(ns)
+    sns = mtkcompile(ns)
     prob = NonlinearProblem(sns, zeros(length(get_unknowns(sns))), [])
 
     sol = NonlinearSolve.solve(prob, NewtonRaphson())