diff --git a/Project.toml b/Project.toml
index 65a904241c..3eb0c183d5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>"]
-version = "3.6.3"
+version = "3.6.4"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 322470c9b3..064129848d 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -63,16 +63,17 @@ function calculate_factorized_W(sys::AbstractODESystem, simplify=true)
     isempty(sys.Wfact[]) || return (sys.Wfact[],sys.Wfact_t[])
 
     jac = calculate_jacobian(sys)
+    M = calculate_massmatrix(sys)
     gam = Variable(:__MTKWgamma)()
 
-    W = - LinearAlgebra.I + gam*jac
+    W = - M + gam*jac
     Wfact = lu(W, Val(false), check=false).factors
 
     if simplify
         Wfact = ModelingToolkit.simplify.(Wfact)
     end
 
-    W_t = - LinearAlgebra.I/gam + jac
+    W_t = - M/gam + jac
     Wfact_t = lu(W_t, Val(false), check=false).factors
     if simplify
         Wfact_t = ModelingToolkit.simplify.(Wfact_t)
diff --git a/test/odesystem.jl b/test/odesystem.jl
index 6813adde69..fcd0e93a94 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -212,9 +212,8 @@ p  = [k₁ => 0.04,
 tspan = (0.0,100000.0)
 prob1 = ODEProblem(sys,u0,tspan,p)
 prob2 = ODEProblem(sys,u0,tspan,p,jac=true)
-# Wfact version is not very stable because of the lack of pivoting
 prob3 = ODEProblem(sys,u0,tspan,p,Wfact=true,Wfact_t=true)
-for (prob, atol) in [(prob1, 1e-12), (prob2, 1e-12), (prob3, 0.1)]
+for (prob, atol) in [(prob1, 1e-12), (prob2, 1e-12), (prob3, 1e-12)]
     sol = solve(prob, Rodas5())
     @test all(x->≈(sum(x), 1.0, atol=atol), sol.u)
 end