diff --git a/Project.toml b/Project.toml
index fefc43e..3093930 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
diff --git a/src/DiffEqsWrappers/DiffEqsWrappers.jl b/src/DiffEqsWrappers/DiffEqsWrappers.jl
new file mode 100644
index 0000000..580a778
--- /dev/null
+++ b/src/DiffEqsWrappers/DiffEqsWrappers.jl
@@ -0,0 +1,91 @@
+module DiffEqWrappers
+
+using Test
+
+using GridapODEs.TransientFETools: TransientFEOperator
+
+using GridapODEs.ODETools: allocate_cache
+using GridapODEs.ODETools: update_cache!
+using GridapODEs.ODETools: residual!
+using GridapODEs.ODETools: jacobian_and_jacobian_t!
+using GridapODEs.ODETools: jacobian!
+using GridapODEs.ODETools: jacobian_t!
+
+using Gridap.Algebra: allocate_jacobian
+
+using Gridap.FESpaces: get_algebraic_operator
+
+using Gridap.Algebra: fill_entries!
+
+export prototype_jacobian
+export prototype_mass
+export prototype_stiffness
+
+export diffeq_wrappers
+
+"""
+  This method takes a `FEOperator` and returns some methods that can be used
+  in `DifferentialEquations.jl`. Assuming we want to solve the nonlinear ODE
+  `res(t,u,du) = 0`, we return:
+  1. `residual!(res, du, u, p, t)`: It returns the residual (in `res`) at
+     `(u,du,t)`, following the signature in `DifferentialEquations`.
+     For the moment, we do not support parameters.
+  2. `jacobian!(jac, du, u, p, gamma, t)`: Idem for the Jacobian. It returns
+     `∂res/∂du*γ + ∂res/∂u`
+  3. `mass!(mass, du, u, p, t)`: Idem for the mass matrix. It returns
+     `∂res/∂du`
+  4. `stiffness!(stif, du, u, p, t)`: Idem for the mass matrix. It returns
+     `∂res/∂u`
+"""
+function diffeq_wrappers(op)
+
+  ode_op = get_algebraic_operator(op)
+  ode_cache = allocate_cache(ode_op)
+
+  function _residual!(res, du, u, p, t)
+    # TO DO (minor): Improve update_cache! st do nothing if same time t as in the cache
+    # now it would be done twice (residual and jacobian)
+    ode_cache = update_cache!(ode_cache, ode_op, t)
+    residual!(res, ode_op, t, u, du, ode_cache)
+  end
+
+  function _jacobian!(jac, du, u, p, gamma, t)
+    ode_cache = update_cache!(ode_cache, ode_op, t)
+    z = zero(eltype(jac))
+    fill_entries!(jac, z)
+    jacobian_and_jacobian_t!(jac, ode_op, t, u, du, gamma, ode_cache)
+  end
+
+  function _mass!(mass, du, u, p, t)
+    ode_cache = update_cache!(ode_cache, ode_op, t)
+    z = zero(eltype(mass))
+    fill_entries!(mass, z)
+    jacobian_t!(mass, ode_op, t, u, du, 1.0, ode_cache)
+  end
+
+  function _stiffness!(stif, du, u, p, t)
+    ode_cache = update_cache!(ode_cache, ode_op, t)
+    z = zero(eltype(stif))
+    fill_entries!(stif, z)
+    jacobian!(stif, ode_op, t, u, du, ode_cache)
+  end
+
+  return _residual!, _jacobian!, _mass!, _stiffness!
+
+end
+
+"""
+  It allocates the Jacobian (or mass or stiffness) matrix, given the `FEOperator`
+  and a vector of size total number of unknowns
+"""
+function prototype_jacobian(op::TransientFEOperator,u0)
+  ode_op = get_algebraic_operator(op)
+  ode_cache = allocate_cache(ode_op) # Not acceptable in terms of performance
+  return allocate_jacobian(ode_op, u0, ode_cache)
+end
+
+const prototype_mass = prototype_jacobian
+
+const prototype_stiffness = prototype_jacobian
+
+end #module
diff --git a/src/GridapODEs.jl b/src/GridapODEs.jl
index df068de..97a096b 100644
--- a/src/GridapODEs.jl
+++ b/src/GridapODEs.jl
@@ -4,6 +4,6 @@ include("ODETools/ODETools.jl")
 
 include("TransientFETools/TransientFETools.jl")
 
-# include("Exports.jl")
+include("DiffEqsWrappers/DiffEqsWrappers.jl")
 
 end #module
diff --git a/test/DiffEqsWrappersTests/DiffEqsTests.jl b/test/DiffEqsWrappersTests/DiffEqsTests.jl
new file mode 100644
index 0000000..36a750f
--- /dev/null
+++ b/test/DiffEqsWrappersTests/DiffEqsTests.jl
@@ -0,0 +1,127 @@
+module DiffEqsWrapperTests
+
+using Test
+using Gridap
+using GridapODEs
+using GridapODEs.ODETools
+using GridapODEs.TransientFETools
+using GridapODEs.DiffEqWrappers
+
+# using DifferentialEquations
+using Sundials
+using Gridap.Algebra: NewtonRaphsonSolver
+using Base.Iterators
+
+# FE problem (heat eq) using Gridap
+function fe_problem(u, n)
+
+  f(t) = x -> ∂t(u)(x, t) - Δ(u(t))(x)
+
+  domain = (0, 1, 0, 1)
+  partition = (n, n)
+  model = CartesianDiscreteModel(domain, partition)
+
+  order = 1
+
+  V0 = FESpace(
+    reffe = :Lagrangian,
+    order = order,
+    valuetype = Float64,
+    conformity = :H1,
+    model = model,
+    dirichlet_tags = "boundary",
+  )
+  U = TransientTrialFESpace(V0, u)
+
+  trian = Triangulation(model)
+  degree = 2 * order
+  quad = CellQuadrature(trian, degree)
+
+  a(u, v) = ∇(v) ⋅ ∇(u)
+  b(v, t) = v * f(t)
+
+  res(t, u, ut, v) = a(u, v) + ut * v - b(v, t)
+  jac(t, u, ut, du, v) = a(du, v)
+  jac_t(t, u, ut, dut, v) = dut * v
+
+  t_Ω = FETerm(res, jac, jac_t, trian, quad)
+  op = TransientFEOperator(U, V0, t_Ω)
+
+  U0 = U(0.0)
+  uh0 = interpolate_everywhere(U0, u(0.0))
+  u0 = get_free_values(uh0)
+
+  return op, u0
+
+end
+
+# Solving the heat equation using GridapODEs and DiffEqs
+tspan = (0.0, 1.0)
+
+u(x, t) = t
+u(t) = x -> u(x, t)
+
+# ISSUE 1: When I choose n > 2, even though the problem that we will solve is
+# linear, the Sundials solvers seems to have convergence issues in the nonlinear
+# solver (?). Ut returns errors
+# [IDAS ERROR]  IDACalcIC Newton/Linesearch algorithm failed to converge.
+# ISSUE 2: When I pass `jac_prototype` the code gets stuck
+
+n = 3 # cells per dim (2D)
+op, u0 = fe_problem(u,n)
+
+# Some checks
+res!, jac!, mass!, stif! = diffeq_wrappers(op)
+J = prototype_jacobian(op,u0)
+r = copy(u0)
+θ = 1.0; t0 = 0.0; tF = 1.0; dt = 0.1; tθ = 1.0; dtθ = dt*θ
+res!(r, u0, u0, nothing, tθ)
+jac!(J, u0, u0, nothing, (1 / dtθ), tθ)
+
+K = prototype_jacobian(op,u0)
+M = prototype_jacobian(op,u0)
+stif!(K, u0, u0, nothing, tθ)
+mass!(M, u0, u0, nothing, tθ)
+# Here you have the mass matrix M
+
+@test (1/dtθ)*M+K ≈ J
+
+# To explore the Sundials solver options, e.g., BE with fixed time step dtd
+f_iip = DAEFunction{true}(res!; jac = jac!)#, jac_prototype=J)
+# jac_prototype is the way to pass my pre-allocated jacobian matrix
+prob_iip = DAEProblem{true}(f_iip, u0, u0, tspan, differential_vars = [true])
+# When I pass `jac_prototype` the code get stuck here:
+sol_iip = Sundials.solve(prob_iip, IDA(), reltol = 1e-8, abstol = 1e-8)
+# @show sol_iip.u
+
+# or iterator version
+integ = init(prob_iip, IDA(), reltol = 1e-8, abstol = 1e-8)
+# step!(integ)
+
+# Show using integrators as iterators
+for i in take(integ, 100)
+  # @show integ.u
+end
+
+end # module
+
+# FUTURE WORK: Check other options, not only Sundials
+
+# ISSUE: Future work, add own (non)linear solver.
+# Let us assume that we just want to consider a fixed point algorithm
+# and we consider an implicit time integration of a nonlinear PDE.
+# Our solvers are efficient since they re-use cache among
+# nonlinear steps (e.g., symbolic factorization, arrays for numerical factorization)
+# and for the transient case in our library, we can also reuse all this
+# between time steps. Could we attain something like this using
+# DifferentialEquations/Sundials?
+# @ChrisRackauckas suggests to take a look at:
+# https://docs.sciml.ai/latest/tutorials/advanced_ode_example/ shows swapping out linear solvers.
+# https://docs.sciml.ai/latest/features/linear_nonlinear/ is all of the extra details.
+
+# Try to pass solver too
+# ls = LUSolver()
+# ls_cache = nothing
+# x = copy(u0)
+# solve!(x,J,r,ls_cache)
+#
diff --git a/test/DiffEqsWrappersTests/runtests.jl b/test/DiffEqsWrappersTests/runtests.jl
new file mode 100644
index 0000000..fd0d0a1
--- /dev/null
+++ b/test/DiffEqsWrappersTests/runtests.jl
@@ -0,0 +1,7 @@
+module DiffEqsWrappersTests
+
+using Test
+
+@testset "DiffEqWrappers" begin include("DiffEqsTests.jl") end
+
+end # module
diff --git a/test/runtests.jl b/test/runtests.jl
index 3a77daf..edc4e51 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,11 +1,13 @@
 module GridapODEsRunTests
 
 using Test
+@time @testset "DiffEqsWrappers" begin include("DiffEqsWrappersTests/runtests.jl") end
 
 @time @testset "ODETools" begin include("ODEsTests/runtests.jl") end
 
 @time @testset "TransientFETools" begin include("TransientFEsTests/runtests.jl") end
 
+
 # include("../bench/runbenchs.jl")
 
 end #module