diff --git a/Project.toml b/Project.toml
index 20f75861..c5c80d2b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,11 +8,14 @@ ACSets = "227ef7b5-1206-438b-ac65-934d6da304b8"
 CombinatorialSpaces = "b1c52339-7909-45ad-8b6a-6e388f7c67f2"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 DiagrammaticEquations = "6f00c28b-6bed-4403-80fa-30e0dc12f317"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
diff --git a/src/operators.jl b/src/operators.jl
index 1aab8c82..70e6194b 100644
--- a/src/operators.jl
+++ b/src/operators.jl
@@ -10,6 +10,7 @@ function default_dec_cu_matrix_generate() end;
 
 function default_dec_matrix_generate(fs::PrimalGeometricMapSeries, my_symbol::Symbol, hodge::DiscreteHodge)
   op = @match my_symbol begin
+    # Inverse Laplacians
     :Δ₀⁻¹ => dec_Δ⁻¹(Val{0}, fs)
     _ => default_dec_matrix_generate(finest_mesh(fs), my_symbol, hodge)
   end
@@ -62,8 +63,11 @@ function default_dec_matrix_generate(sd::HasDeltaSet, my_symbol::Symbol, hodge::
     :ℒ₁ => ℒ_dd(Tuple{1,1}, sd)
 
     # Dual Laplacians
-    :Δᵈ₀ => Δᵈ(Val{0},sd)
-    :Δᵈ₁ => Δᵈ(Val{1},sd)
+    :Δᵈ₀ => Δᵈ(Val{0}, sd)
+    :Δᵈ₁ => Δᵈ(Val{1}, sd)
+
+    # Inverse Laplacians
+    :Δ₀⁻¹ => dec_inv_lap_solver(Val{0}, sd)
 
     # Musical Isomorphisms
     :♯ => dec_♯_p(sd)
@@ -154,6 +158,11 @@ function dec_avg₀₁(sd::HasDeltaSet)
   (avg_mat, x -> avg_mat * x)
 end
 
+function dec_inv_lap_solver(::Type{Val{0}}, sd::HasDeltaSet)
+  inv_lap = LinearAlgebra.factorize(∇²(0, sd))
+  x -> inv_lap \ x
+end
+
 function default_dec_generate(sd::HasDeltaSet, my_symbol::Symbol, hodge::DiscreteHodge=GeometricHodge())
 
   op = @match my_symbol begin
diff --git a/test/simulation.jl b/test/simulation.jl
index 01a27437..a69f2521 100644
--- a/test/simulation.jl
+++ b/test/simulation.jl
@@ -770,13 +770,21 @@ end
     end
   end
 
-  sim = eval(gensim(inv_lap; multigrid=true))
+  sim = eval(gensim(inv_lap))
+  sim_mg = eval(gensim(inv_lap; multigrid=true))
+
+  f = sim(our_mesh, generate);
+  f_mg = sim_mg(series, generate);
 
-  f = sim(series, generate);
   u = ComponentArray(U=b)
   du = similar(u)
 
+  # Regular mesh
   f(du, u, 0, ())
+  @test norm(lap*du.U-b)/norm(b) < 1e-15
+
+  # Multigrid
+  f_mg(du, u, 0, ())
   @test norm(lap*du.U-b)/norm(b) < 1e-6
 end