diff --git a/Project.toml b/Project.toml
index d32e1742a..db8818b1b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -26,13 +26,13 @@ DomainSets = "0.5, 0.6"
 IfElse = "0.1"
 Interpolations = "0.14"
 Latexify = "0.15"
-ModelingToolkit = "8.22"
+ModelingToolkit = "8"
 OrdinaryDiffEq = "6"
 SafeTestsets = "0.0.1"
 SciMLBase = "1.57"
 StaticArrays = "1"
-SymbolicUtils = "0.19, 1"
-Symbolics = "4, 5"
+SymbolicUtils = "1"
+Symbolics = "5"
 TermInterface = "0.2, 0.3"
 julia = "1.6"
 
diff --git a/docs/src/tutorials/PIDE.md b/docs/src/tutorials/PIDE.md
index a553fe1bd..d7ff40226 100644
--- a/docs/src/tutorials/PIDE.md
+++ b/docs/src/tutorials/PIDE.md
@@ -56,6 +56,8 @@ To have an integral over the whole domain, be sure to wrap the integral in an au
 Due to a limitation, the whole domain integral needs to have the same arguments as the integrand, but is constant in x. To use it in an equation one dimension lower, use a boundary value like integral(t, 0.0)
 
 ```@example integrals2
+using MethodOfLines, ModelingToolkit, DomainSets
+
 @parameters t, x
 @variables integrand(..) integral(..)
 Dt = Differential(t)
diff --git a/src/MOL_symbolic_utils.jl b/src/MOL_symbolic_utils.jl
index 4695cddaa..74b2913f1 100644
--- a/src/MOL_symbolic_utils.jl
+++ b/src/MOL_symbolic_utils.jl
@@ -90,17 +90,13 @@ subs_alleqs!(eqs, rules) = map!(eq -> substitute(eq.lhs, rules) ~ substitute(eq.
 find all the dependent variables given by depvar_ops in an expression
 """
 function get_depvars(eq, depvar_ops)
-    S = Symbolics
-    SU = SymbolicUtils
     depvars = Set()
-    if eq isa Num
-        eq = eq.val
-    end
-    if S.istree(eq)
-        if eq isa Term && any(u -> isequal(operation(eq), u), depvar_ops)
+    eq = safe_unwrap(eq)
+    if istree(eq)
+        if any(u -> isequal(operation(eq), u), depvar_ops)
             push!(depvars, eq)
         else
-            for o in map(x -> get_depvars(x, depvar_ops), SU.arguments(eq))
+            for o in map(x -> get_depvars(x, depvar_ops), arguments(eq))
                 union!(depvars, o)
             end
         end
diff --git a/src/discretization/discretize_vars.jl b/src/discretization/discretize_vars.jl
index 9e57a6b19..6f75312d3 100644
--- a/src/discretization/discretize_vars.jl
+++ b/src/discretization/discretize_vars.jl
@@ -129,7 +129,7 @@ function DiscreteSpace(vars, discretization::MOLFiniteDifference{G}) where {G}
 
     depvarsdisc = map(depvars) do u
         op = SymbolicUtils.operation(u)
-        if op isa SymbolicUtils.Term{SymbolicUtils.FnType{Tuple,Real},Nothing}
+        if op isa  SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}
             sym = Symbol(string(op))
         else
             sym = nameof(op)
diff --git a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl
index 32fc91a60..d97ca7886 100644
--- a/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl
+++ b/src/discretization/schemes/nonlinear_laplacian/nonlinear_laplacian.jl
@@ -72,13 +72,13 @@ function cartesian_nonlinear_laplacian(expr, II, derivweights, s::DiscreteSpace,
 end
 
 @inline function generate_nonlinlap_rules(II::CartesianIndex, s::DiscreteSpace, depvars, derivweights::DifferentialDiscretizer, bcmap, indexmap, terms)
-    rules = reduce(safe_vcat, [vec([@rule *(~~c, $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)), ~~d) => *(~c..., cartesian_nonlinear_laplacian(*(a..., b...), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u), ~d...) for x in params(u, s)]) for u in depvars], init = [])
+    rules = reduce(safe_vcat, [vec([@rule *(~~c, $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)), ~~d) => *(~c..., cartesian_nonlinear_laplacian(*(~a..., ~b...), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u), ~d...) for x in params(u, s)]) for u in depvars], init = [])
 
-    rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)) => cartesian_nonlinear_laplacian(*(a..., b...), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u) for x in params(u, s)]) for u in depvars], init = []))
+    rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule $(Differential(x))(*(~~a, $(Differential(x))(u), ~~b)) => cartesian_nonlinear_laplacian(*(~a..., ~b...), Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u) for x in params(u, s)]) for u in depvars], init = []))
 
     rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule ($(Differential(x))($(Differential(x))(u) / ~a)) => cartesian_nonlinear_laplacian(1 / ~a, Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u) for x in params(u, s)]) for u in depvars], init = []))
 
-    rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule *(~~b, ($(Differential(x))($(Differential(x))(u) / ~a)), ~~c) => *(b..., c..., cartesian_nonlinear_laplacian(1 / ~a, Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u)) for x in params(u, s)]) for u in depvars], init = []))
+    rules = safe_vcat(rules, reduce(safe_vcat, [vec([@rule *(~~b, ($(Differential(x))($(Differential(x))(u) / ~a)), ~~c) => *(~b..., ~c..., cartesian_nonlinear_laplacian(1 / ~a, Idx(II, s, u, indexmap), derivweights, s, filter_interfaces(bcmap[operation(u)][x]), depvars, x, u)) for x in params(u, s)]) for u in depvars], init = []))
 
     nonlinlap_rules = []
     for t in terms
diff --git a/test/components/interiormap_test.jl b/test/components/interiormap_test.jl
index c6aa11673..57a9a84c3 100644
--- a/test/components/interiormap_test.jl
+++ b/test/components/interiormap_test.jl
@@ -38,9 +38,9 @@ const bigint = div(typemax(Int), 2)
 
     m = MethodOfLines.buildmatrix(pde, s)
     if VERSION >= v"1.7"
-        @test m == [2+bigint 0 1; 0 2+bigint 1; 1 1 2+bigint] # Test the matrix is the identity matrix
-    else
         @test m == [1 2+bigint 0; 1 0 2+bigint; 2+bigint 1 1]
+    else
+        @test m == [0 2+bigint 1; 2+bigint 0 1; 1 1 2+bigint] # Test the matrix is the identity matrix
     end
 
 end
@@ -78,9 +78,9 @@ end
     s = MethodOfLines.DiscreteSpace(v, disc)
     m = MethodOfLines.buildmatrix(pde, s)
     if VERSION >= v"1.7"
-        @test m == [2 0 2; 0 3 3; 4 4 4] # Test the matrix is the identity matrix
-    else
         @test m == [2 2 0; 3 0 3; 4 4 4]
+    else
+        @test m == [0 2 2; 3 0 3; 4 4 4]
     end
 end
 #
@@ -117,9 +117,9 @@ end
 
     m = MethodOfLines.buildmatrix(pde, s)
     if VERSION >= v"1.7"
-        @test m == [2 0 1; 0 3 3; 5 5 4] # Test the matrix is the identity matrix
-    else
         @test m == [1 2 0; 3 0 3; 4 5 5]
+    else
+        @test m == [0 2 1; 3 0 3; 5 5 4]
     end
 
 end
@@ -157,9 +157,9 @@ end
     s = MethodOfLines.DiscreteSpace(v, disc)
     m = MethodOfLines.buildmatrix(pde, s)
     if VERSION >= v"1.7"
-        @test m == [2 2 0; 1 0 1; 0 1 2]
-    else
         @test m == [0 2 2; 1 1 0; 2 0 1]
+    else
+        @test m == [2 2 0; 0 1 1; 1 0 2]
     end
 
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 1b7597110..05d2aa01b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,6 +7,25 @@ const is_TRAVIS = haskey(ENV, "TRAVIS")
 # Start Test Script
 
 @time begin
+    if GROUP == "All" || GROUP == "Components"
+        @time @safetestset "MOLFiniteDifference Utils" begin
+            include("utils_test.jl")
+        end
+        @time @safetestset "Discretization of space and grid types" begin
+            include("components/DiscreteSpace.jl")
+        end
+        @time @safetestset "Variable PDE mapping and interior construction" begin
+            include("components/interiormap_test.jl")
+        end
+        @time @safetestset "Fornberg" begin
+            include("components/MOLfornberg_weights.jl")
+        end
+        @time @safetestset "ODEFunction" begin
+            include("components/ODEFunction_test.jl")
+        end
+        #@time @safetestset "Finite Difference Schemes" begin include("components/finite_diff_schemes.jl") end
+    end
+
     if GROUP == "All" || GROUP == "Integrals"
         @time @safetestset "MOLFiniteDifference Interface: Integrals" begin
             include("pde_systems/MOL_1D_Integration.jl")
@@ -31,25 +50,6 @@ const is_TRAVIS = haskey(ENV, "TRAVIS")
         end
     end
 
-    if GROUP == "All" || GROUP == "Components"
-        #@time @safetestset "Test for regression against original code" begin include("regression_test.jl") end
-        @time @safetestset "MOLFiniteDifference Utils" begin
-            include("utils_test.jl")
-        end
-        @time @safetestset "Discretization of space and grid types" begin
-            include("components/DiscreteSpace.jl")
-        end
-        @time @safetestset "Variable PDE mapping and interior construction" begin
-            include("components/interiormap_test.jl")
-        end
-        @time @safetestset "Fornberg" begin
-            include("components/MOLfornberg_weights.jl")
-        end
-        @time @safetestset "ODEFunction" begin
-            include("components/ODEFunction_test.jl")
-        end
-        #@time @safetestset "Finite Difference Schemes" begin include("components/finite_diff_schemes.jl") end
-    end
 
     if GROUP == "All" || GROUP == "Nonlinear_Diffusion"
         @time @safetestset "MOLFiniteDifference Interface: 1D Non-Linear Diffusion" begin