diff --git a/src/GridapODEs/ODETools/AffineNewmark.jl b/src/GridapODEs/ODETools/AffineNewmark.jl
index 905af110d..7e229e71c 100644
--- a/src/GridapODEs/ODETools/AffineNewmark.jl
+++ b/src/GridapODEs/ODETools/AffineNewmark.jl
@@ -97,7 +97,7 @@ function jacobian!(A::AbstractMatrix,op::NewmarkAffineOperator,x::AbstractVector
   a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0
   v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0
   z = zero(eltype(A))
-  fill!(A,z)
+  fillstored!(A,z)
   jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache)
 end
 
diff --git a/src/GridapODEs/ODETools/AffineThetaMethod.jl b/src/GridapODEs/ODETools/AffineThetaMethod.jl
index afa110ce6..15ae7e7f7 100644
--- a/src/GridapODEs/ODETools/AffineThetaMethod.jl
+++ b/src/GridapODEs/ODETools/AffineThetaMethod.jl
@@ -106,13 +106,13 @@ end
 
 function _matrix!(A,odeop,tθ,dtθ,u0,ode_cache,vθ)
   z = zero(eltype(A))
-  fill!(A,z)
+  fillstored!(A,z)
   jacobians!(A,odeop,tθ,(vθ,vθ),(1.0,1/dtθ),ode_cache)
 end
 
 function _mass_matrix!(A,odeop,tθ,dtθ,u0,ode_cache,vθ)
   z = zero(eltype(A))
-  fill!(A,z)
+  fillstored!(A,z)
   jacobian!(A,odeop,tθ,(vθ,vθ),2,(1/dtθ),ode_cache)
 end
 
@@ -143,7 +143,7 @@ function ThetaMethodConstantOperator(odeop::ConstantODEOperator,tθ::Float64,dt
   residual!(b,odeop,tθ,(u0,vθ),ode_cache)
   b = -1*b
   z = zero(eltype(A))
-  fill!(A,z)
+  fillstored!(A,z)
   jacobians!(A,odeop,tθ,(vθ,vθ),(1.0,1/dtθ),ode_cache)
   return A, b
 end
diff --git a/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl b/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl
index d40a55bc9..06699bcf0 100644
--- a/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl
+++ b/src/GridapODEs/ODETools/ConstantMatrixNewmark.jl
@@ -92,7 +92,7 @@ function jacobian!(A::AbstractMatrix,op::NewmarkConstantMatrixOperator,x::Abstra
   a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0
   v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0
   z = zero(eltype(A))
-  fill!(A,z)
+  fillstored!(A,z)
   jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache)
 end
 
diff --git a/src/GridapODEs/ODETools/ConstantNewmark.jl b/src/GridapODEs/ODETools/ConstantNewmark.jl
index 1b32f67a8..685a117df 100644
--- a/src/GridapODEs/ODETools/ConstantNewmark.jl
+++ b/src/GridapODEs/ODETools/ConstantNewmark.jl
@@ -125,7 +125,7 @@ function jacobian!(A::AbstractMatrix,op::NewmarkConstantOperator,x::AbstractVect
   a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0
   v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0
   z = zero(eltype(A))
-  fill!(A,z)
+  fillstored!(A,z)
   jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache)
 end
 
@@ -136,7 +136,7 @@ function _mass_matrix!(A::AbstractMatrix,op::NewmarkConstantOperator,x::Abstract
   a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0
   v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0
   z = zero(eltype(A))
-  fill!(A,z)
+  fillstored!(A,z)
   jacobian!(A,op.odeop,op.t1,(u1,v1,a1),3,1.0,cache)
 end
 
@@ -147,6 +147,6 @@ function _damping_matrix!(A::AbstractMatrix,op::NewmarkConstantOperator,x::Abstr
   a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0
   v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0
   z = zero(eltype(A))
-  fill!(A,z)
+  fillstored!(A,z)
   jacobian!(A,op.odeop,op.t1,(u1,v1,a1),2,1.0,cache)
 end
diff --git a/src/GridapODEs/ODETools/ForwardEuler.jl b/src/GridapODEs/ODETools/ForwardEuler.jl
index 7b68b2f72..b177b55b6 100644
--- a/src/GridapODEs/ODETools/ForwardEuler.jl
+++ b/src/GridapODEs/ODETools/ForwardEuler.jl
@@ -59,7 +59,7 @@ function jacobian!(A::AbstractMatrix,op::ForwardEulerNonlinearOperator,x::Abstra
   vf = op.vf
   vf = (x-op.u0)/op.dt
   z = zero(eltype(A))
-  fill!(A,z)
+  fillstored!(A,z)
   jacobians!(A,op.odeop,op.tf,(op.u0,vf),(0,1/op.dt),op.ode_cache)
 end
 
diff --git a/src/GridapODEs/ODETools/Newmark.jl b/src/GridapODEs/ODETools/Newmark.jl
index 8588c0ce4..eda2fe93f 100644
--- a/src/GridapODEs/ODETools/Newmark.jl
+++ b/src/GridapODEs/ODETools/Newmark.jl
@@ -75,7 +75,7 @@ function jacobian!(A::AbstractMatrix,op::NewmarkNonlinearOperator,x::AbstractVec
   a1 = 1.0/(op.β*op.dt^2)*(u1-u0) - 1.0/(op.β*op.dt)*v0 - (1-2*op.β)/(2*op.β)*a0
   v1 = op.γ/(op.β*op.dt)*(u1-u0) + (1-op.γ/op.β)*v0 + op.dt*(1-op.γ/(2*op.β))*a0
   z = zero(eltype(A))
-  fill!(A,z)
+  fillstored!(A,z)
   jacobians!(A,op.odeop,op.t1,(u1,v1,a1),(1.0,op.γ/(op.β*op.dt),1.0/(op.β*op.dt^2)),cache)
 end
 
diff --git a/src/GridapODEs/ODETools/ODETools.jl b/src/GridapODEs/ODETools/ODETools.jl
index 9e4924796..bb06d8b1a 100644
--- a/src/GridapODEs/ODETools/ODETools.jl
+++ b/src/GridapODEs/ODETools/ODETools.jl
@@ -10,7 +10,8 @@ using Test
 using DocStringExtensions
 
 using ForwardDiff
-using LinearAlgebra: fill!
+using LinearAlgebra: fillstored!
+using SparseArrays: issparse
 
 const ϵ = 100*eps()
 export ∂t
diff --git a/src/GridapODEs/ODETools/RungeKutta.jl b/src/GridapODEs/ODETools/RungeKutta.jl
index 53a2ad900..ee2768380 100644
--- a/src/GridapODEs/ODETools/RungeKutta.jl
+++ b/src/GridapODEs/ODETools/RungeKutta.jl
@@ -175,7 +175,7 @@ function jacobian!(A::AbstractMatrix,op::RungeKuttaNonlinearOperator,x::Abstract
   vi = op.vi
   vi = (x-op.u0)/(op.a[op.i,op.i]*op.dt)
   z = zero(eltype(A))
-  fill!(A,z)
+  fillstored!(A,z)
   jacobians!(A,op.odeop,op.ti,(ui,vi),(1.0,1.0/(op.a[op.i,op.i]*op.dt)),op.ode_cache)
 end
 
diff --git a/src/GridapODEs/ODETools/ThetaMethod.jl b/src/GridapODEs/ODETools/ThetaMethod.jl
index 2bc772d53..dae23eb32 100644
--- a/src/GridapODEs/ODETools/ThetaMethod.jl
+++ b/src/GridapODEs/ODETools/ThetaMethod.jl
@@ -79,7 +79,7 @@ function jacobian!(A::AbstractMatrix,op::ThetaMethodNonlinearOperator,x::Abstrac
   vθ = op.vθ
   vθ = (x-op.u0)/op.dtθ
   z = zero(eltype(A))
-  fill!(A,z)
+  fillstored!(A,z)
   jacobians!(A,op.odeop,op.tθ,(uF,vθ),(1.0,1/op.dtθ),op.ode_cache)
 end
 
diff --git a/test/GridapODEsTests/ODEsTests/ODEOperatorMocks.jl b/test/GridapODEsTests/ODEsTests/ODEOperatorMocks.jl
index ba7453901..5635aa7c6 100644
--- a/test/GridapODEsTests/ODEsTests/ODEOperatorMocks.jl
+++ b/test/GridapODEsTests/ODEsTests/ODEOperatorMocks.jl
@@ -16,6 +16,7 @@ import Gridap.GridapODEs.ODETools: jacobian!
 import Gridap.GridapODEs.ODETools: jacobians!
 import Gridap.GridapODEs.ODETools: allocate_jacobian
 import Gridap.GridapODEs.ODETools: residual!
+using SparseArrays: spzeros
 
 struct ODEOperatorMock{T<:Real,C} <: ODEOperator{C}
   a::T
@@ -104,7 +105,7 @@ function jacobians!(
 end
 
 function allocate_jacobian(op::ODEOperatorMock,u::AbstractVector,cache)
-  zeros(2,2)
+  spzeros(2,2)
 end
 
 allocate_cache(op::ODEOperatorMock) = nothing
diff --git a/test/GridapODEsTests/ODEsTests/ODESolverMocks.jl b/test/GridapODEsTests/ODEsTests/ODESolverMocks.jl
index b65c5dc12..7ae341376 100644
--- a/test/GridapODEsTests/ODEsTests/ODESolverMocks.jl
+++ b/test/GridapODEsTests/ODEsTests/ODESolverMocks.jl
@@ -60,7 +60,7 @@ end
 function solve!(x::AbstractVector,nls::NLSolverMock,nlop::NonlinearOperator,cache::Nothing)
   r = residual(nlop,x)
   J = jacobian(nlop,x)
-  dx = inv(J)*(-r)
+  dx = inv(Matrix(J))*(-r)
   x.= x.+dx
   cache = (r,J,dx)
 end
@@ -69,7 +69,7 @@ function solve!(x::AbstractVector,nls::NLSolverMock,nlop::NonlinearOperator,cach
   r, J, dx = cache
   residual!(r, nlop, x)
   jacobian!(J, nlop, x)
-  dx = inv(J)*(-r)
+  dx = inv(Matrix(J))*(-r)
   x.= x.+dx
 end