diff --git a/.gitignore b/.gitignore
index 2db1a2f8d..23e1064ba 100644
--- a/.gitignore
+++ b/.gitignore
@@ -12,7 +12,7 @@ test/plots
 test/plots/*
 dev
 dev/*
-
+*.gif
 # Build artifacts for creating documentation generated by the Documenter package
 docs/build/
 docs/site/
diff --git a/src/system_parsing/interior_map.jl b/src/system_parsing/interior_map.jl
index 2c7e9f243..fe0cad034 100644
--- a/src/system_parsing/interior_map.jl
+++ b/src/system_parsing/interior_map.jl
@@ -31,8 +31,8 @@ function PDEBase.construct_var_equation_mapping(
         pdes::Vector{Equation}, boundarymap, s::DiscreteSpace{N, M},
         discretization::MOLFiniteDifference) where {N, M}
     @assert length(pdes)==M "There must be the same number of equations and unknowns, got $(length(pdes)) equations and $(M) unknowns"
-    m = buildmatrix(pdes, s)
-    varmap = Dict(build_variable_mapping(m, s.ū, pdes))
+    m, var_dims, varpalsmap = buildmatrix(pdes, s)
+    varmap = Dict(build_variable_mapping(m, s.ū, pdes))              
 
     # Determine the interiors for each pde
     vlower = []
@@ -41,32 +41,41 @@ function PDEBase.construct_var_equation_mapping(
 
     interior = map(pdes) do pde
         u = varmap[pde]
-        boundaries = mapreduce(x -> boundarymap[operation(u)][x], vcat, s.x̄)
-        n = ndims(u, s)
-        lower = zeros(Int, n)
-        upper = zeros(Int, n)
-        # Determine thec number of points to remove from each end of the domain for each dimension
-        for b in boundaries
-            #@show b
-            clip_interior!!(lower, upper, s, b)
-        end
-        push!(vlower, pde => lower)
-        push!(vupper, pde => upper)
-        #TODO: Allow assymmetry
-        pdeorders = Dict(map(x -> x => d_orders(x, [pde]), s.x̄))
-
-        # Add ghost points to pad stencil extents
-        lowerextents, upperextents = calculate_stencil_extents(
-            s, u, discretization, pdeorders, boundarymap)
-        push!(extents, pde => (lowerextents, upperextents))
-        lower = [max(e, l) for (e, l) in zip(lowerextents, lower)]
-        upper = [max(e, u) for (e, u) in zip(upperextents, upper)]
-        mindomsize = lower .+ upper .+ 1
-        if any(tup -> mindomsize[tup[1]] > length(s, tup[2]), enumerate(ivs(u, s)))
-            error("The domain is too small to support the requested discretization, got domain size of $(size(s)).")
+        f = v -> findmax(map(p -> p[2], varpalsmap[v]))[2]
+        if var_dims[u] == 0 && !isnothing(f(u))
+            i = f(u)
+            u = first(varpalsmap[u][i])
+            args = remove(arguments(u), s.time)
+
+            pde => s.Igrid[u][[1:length(s.grid[x]) for x in args]...]
+        else
+            boundaries = mapreduce(x -> boundarymap[operation(u)][x], vcat, s.x̄)
+            n = ndims(u, s)
+            lower = zeros(Int, n)
+            upper = zeros(Int, n)
+            # Determine thec number of points to remove from each end of the domain for each dimension
+            for b in boundaries
+                #@show b
+                clip_interior!!(lower, upper, s, b)
+            end
+            push!(vlower, pde => lower)
+            push!(vupper, pde => upper)
+            #TODO: Allow assymmetry
+            pdeorders = Dict(map(x -> x => d_orders(x, [pde]), s.x̄))
+
+            # Add ghost points to pad stencil extents
+            lowerextents, upperextents = calculate_stencil_extents(
+                s, u, discretization, pdeorders, boundarymap)
+            push!(extents, pde => (lowerextents, upperextents))
+            lower = [max(e, l) for (e, l) in zip(lowerextents, lower)]
+            upper = [max(e, u) for (e, u) in zip(upperextents, upper)]
+            mindomsize = lower .+ upper .+ 1
+            if any(tup -> mindomsize[tup[1]] > length(s, tup[2]), enumerate(ivs(u, s)))
+                error("The domain is too small to support the requested discretization, got domain size of $(size(s)).")
+            end
+            # Don't update this x2i, it is correct.
+            pde => generate_interior(lower, upper, u, s, discretization)
         end
-        # Don't update this x2i, it is correct.
-        pde => generate_interior(lower, upper, u, s, discretization)
     end
 
     pdemap = [k.second => k.first for k in varmap]
@@ -117,26 +126,43 @@ end
 function buildmatrix(pdes, s::DiscreteSpace{N, M}) where {N, M}
     m = zeros(Int, M, M)
     elegiblevars = [getvarmap(pde, s) for pde in pdes]
-    u2i = Dict([u => k for (k, u) in enumerate(s.ū)])
-    #@show elegiblevars, s.ū
+    u2i = Dict([u => k for (k, u) in enumerate(s.ū)])
+    var_dims = Dict([var => variable_dimensionality(var, s) for var in s.ū])
+    
     for (i, varmap) in enumerate(elegiblevars)
-        for var in keys(varmap)
-            m[i, u2i[var]] = varmap[var]
+        for var in keys(varmap)            
+            m[i, u2i[var]] = var_dims[var] == 0 ? varmap[var] - div(typemax(Int), 2) : varmap[var]
         end
     end
-    return m
+    varpalmap = map(filter(x -> var_dims[x] == 0, s.ū)) do u
+        #find i where var_dims[u] shares a row with a non-zero var_dims[v]
+        varpals = []
+        for j in 1:M
+            if m[j, u2i[u]] != 0
+                for k in setdiff(1:M, [u2i[u]])
+                    if m[j, k] != 0 && var_dims[s.ū[k]] > 0
+                        push!(varpals, s.ū[k] => m[j, k])
+                    end
+                end
+            end
+        end
+        u => varpals
+    end |> Dict
+    return m, var_dims, varpalmap
 end
 
-function build_variable_mapping(m, vars, pdes)
+function build_variable_mapping(m::Matrix{Int}, vars::Vector, pdes::Vector)
     notzero(x) = x > 0 ? 1 : 0
     varpdemap = []
     N = length(pdes)
+    
     rows = sum(m, dims = 2)
     cols = sum(m, dims = 1)
     i = findfirst(isequal(0), rows)
     j = findfirst(isequal(0), cols)
     @assert i===nothing "Equation $(pdes[i[1]]) is not an equation for any of the dependent variables."
     @assert j===nothing "Variable $(vars[j[2]]) does not appear in any equation, therefore cannot be solved for"
+    
     for k in 1:N
         # Check if any of the pdes only have one valid variable
         m_ones = notzero.(m)
@@ -240,3 +266,8 @@ function getvarmap(pde, s)
     end
     return varmap
 end
+
+function variable_dimensionality(var, s::DiscreteSpace)
+    args = arguments(var)
+    return count(arg -> any(isequal(arg), s.x̄), args)
+end
diff --git a/test/components/interiormap_test.jl b/test/components/interiormap_test.jl
index a080642ae..60230df67 100644
--- a/test/components/interiormap_test.jl
+++ b/test/components/interiormap_test.jl
@@ -36,7 +36,7 @@ const bigint = div(typemax(Int), 2)
 
     s = MethodOfLines.construct_discrete_space(v, disc)
 
-    m = MethodOfLines.buildmatrix(pde, s)
+    m, _, _ = MethodOfLines.buildmatrix(pde, s)
     
     test = [1 2+bigint 0; 1 0 2+bigint; 2+bigint 1 1]
     perms = permutations([1, 2, 3])
@@ -77,7 +77,7 @@ end
     v = MethodOfLines.VariableMap(pdesys, disc)
 
     s = MethodOfLines.construct_discrete_space(v, disc)
-    m = MethodOfLines.buildmatrix(pde, s)
+    m, _, __ = MethodOfLines.buildmatrix(pde, s)
     test = [2 2 0; 3 0 3; 4 4 4]
     
     perms = permutations([1, 2, 3])
@@ -117,7 +117,7 @@ end
     v = MethodOfLines.VariableMap(pdesys, disc)
     s = MethodOfLines.construct_discrete_space(v, disc)
 
-    m = MethodOfLines.buildmatrix(pde, s)
+    m, _, __ = MethodOfLines.buildmatrix(pde, s)
     test = [1 2 0; 3 0 3; 4 5 5]
     perms = permutations([1, 2, 3])
     @test any(perms) do perm
@@ -157,7 +157,7 @@ end
     v = MethodOfLines.VariableMap(pdesys, disc)
 
     s = MethodOfLines.construct_discrete_space(v, disc)
-    m = MethodOfLines.buildmatrix(pde, s)
+    m, _, __ = MethodOfLines.buildmatrix(pde, s)
     test = [0 2 2; 1 1 0; 2 0 1]
     
     perms = permutations([1, 2, 3])
@@ -167,6 +167,50 @@ end
 
 end
 
+@testset "Test 00d: ignore time-only variables when higher-dimensional variables exist" begin
+    @parameters x, t
+    @variables u(..), v(..), w(..)
+
+    Dxx = Differential(x)^2
+    Dt = Differential(t)
+
+    t_min = 0.
+    t_max = 2.0
+    x_min = 0.
+    x_max = 20.0
+
+    dx = 1.0
+
+    domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)]
+
+    # u and w are functions of t and x, v is a function of t only
+    pde = [Dt(u(t,x)) ~ Dxx(u(t,x)) +  w(t,x),
+           Dt(v(t)) ~ u(t,x) + w(t,x),
+           Dt(w(t,x)) ~ Dxx(w(t,x)) + u(t,x)]
+
+    bcs = [u(0,x) ~ 0, u(t,0) ~ 0, u(t,x_max) ~ 0,
+           w(0,x) ~ 0, w(t,0) ~ 0, w(t,x_max) ~ 0]
+
+    @named pdesys = PDESystem(pde, bcs, domains, [t,x], [u(t,x), v(t), w(t,x)])
+
+    # Test centered order
+    disc = MOLFiniteDifference([x=>dx], t)
+
+    v = MethodOfLines.VariableMap(pdesys, disc)
+    s = MethodOfLines.construct_discrete_space(v, disc)
+
+    m, _, __ = MethodOfLines.buildmatrix(pde, s)
+    
+    # The expected result should ignore v(t) in the mapping process
+    expected = [-4611686018427387903 1 4611686018427387905; 2 1 1; -4611686018427387903 4611686018427387905 1]
+    
+    @test m == expected
+
+    varmap = Dict(MethodOfLines.build_variable_mapping(m, s.ū, pde))
+    
+    @test_broken prob = discretize(pdesys, disc)
+end
+
 @testset "Test 01a: Build variable mapping - one right choice simple" begin
     m = hcat([0, 1, 0],
              [0, 0, 1],
@@ -241,3 +285,4 @@ end
         @test e isa AssertionError
     end
 end
+
diff --git a/test/pde_systems/nonlinear_laplacian_advanced.jl b/test/pde_systems/nonlinear_laplacian_advanced.jl
index c8290b5e2..6ec8aec70 100644
--- a/test/pde_systems/nonlinear_laplacian_advanced.jl
+++ b/test/pde_systems/nonlinear_laplacian_advanced.jl
@@ -32,7 +32,8 @@ function asf(dt, dx, dy)
     halfar_dome(dt, dx, dy, R0, H0, 917)
 end
 
-@test_broken begin#@testset "Halfar ice dome glacier model." begin
+
+@testset "Halfar ice dome glacier model." begin
     rmax = 2 * 1000
     rmin = -rmax
 
@@ -78,11 +79,19 @@ end
 
     sol = solve(prob, FBDF())
 
-    @test SciMLBase.successful_retcode(sol)
+    @test_broken SciMLBase.successful_retcode(sol)
+
 
     solx = sol[x]
     soly = sol[y]
     solt = sol[t]
+    solH = sol[H(t, x, y)]
+
+    #anim = @animate for (i, t) in enumerate(t)
+    #    heatmap(solx, soly, solH[i, :, :])
+    #end
+    #gif(anim, "halfar.gif", fps = 10)
+    
 
     solexact = [asf(unwrap(dt), unwrap(dx), unwrap(dy))
                 for dt in solt, dx in solx, dy in soly]