diff --git a/results.vtu b/results.vtu
deleted file mode 100644
index 6a70a18..0000000
Binary files a/results.vtu and /dev/null differ
diff --git a/src/ODETools/ForwardEuler.jl b/src/ODETools/ForwardEuler.jl
index 67b67da..b177b55 100644
--- a/src/ODETools/ForwardEuler.jl
+++ b/src/ODETools/ForwardEuler.jl
@@ -56,12 +56,11 @@ function residual!(b::AbstractVector,op::ForwardEulerNonlinearOperator,x::Abstra
 end
 
 function jacobian!(A::AbstractMatrix,op::ForwardEulerNonlinearOperator,x::AbstractVector)
-  uF = x
   vf = op.vf
   vf = (x-op.u0)/op.dt
   z = zero(eltype(A))
   fillstored!(A,z)
-  jacobian!(A,op.odeop,op.tf,(op.u0,vf),1,(1/op.dt),op.ode_cache)
+  jacobians!(A,op.odeop,op.tf,(op.u0,vf),(0,1/op.dt),op.ode_cache)
 end
 
 function allocate_residual(op::ForwardEulerNonlinearOperator,x::AbstractVector)
diff --git a/src/TransientFETools/TransientFEOperators.jl b/src/TransientFETools/TransientFEOperators.jl
index 02b4f7f..6705884 100644
--- a/src/TransientFETools/TransientFEOperators.jl
+++ b/src/TransientFETools/TransientFEOperators.jl
@@ -289,7 +289,9 @@ function jacobians!(
   cache)
   _matdata = ()
   for i in 1:get_order(op)+1
-    _matdata = (_matdata...,matdata_jacobian(op,t,xh,i,γ[i]))
+    if (γ[i] > 0.0)
+      _matdata = (_matdata...,matdata_jacobian(op,t,xh,i,γ[i]))
+    end
   end
   matdata = vcat_matdata(_matdata)
   assemble_matrix_add!(A,op.assem_t, matdata)
diff --git a/test/TransientFEsTests/ForwardEulerHeatEquationTests.jl b/test/TransientFEsTests/ForwardEulerHeatEquationTests.jl
index 6cf75a3..759f014 100644
--- a/test/TransientFEsTests/ForwardEulerHeatEquationTests.jl
+++ b/test/TransientFEsTests/ForwardEulerHeatEquationTests.jl
@@ -48,10 +48,9 @@ a(u,v) = ∇(v)⋅∇(u)
 b(v,t) = v*f(t)
 
 res(t,(u,ut),v) = ∫( a(u,v) + ut*v - b(v,t) )dΩ
-jac(t,(u,ut),du,v) = a(du,v)
+jac(t,(u,ut),du,v) = ∫( a(du,v) )dΩ
 jac_t(t,(u,ut),dut,v) = ∫( dut*v )dΩ
 
-#t_Ω = FETerm(res,jac,jac_t,trian,quad)
 op = TransientFEOperator(res,jac,jac_t,U,V0)
 
 t0 = 0.0
@@ -62,8 +61,6 @@ U0 = U(0.0)
 uh0 = interpolate_everywhere(u(0.0),U0)
 
 ls = LUSolver()
-using Gridap.Algebra: NewtonRaphsonSolver
-nls = NLSolver(ls;show_trace=true,method=:newton) #linesearch=BackTracking())
 ode_solver = ThetaMethod(ls,dt,θ)
 
 sol_t = solve(ode_solver,op,uh0,t0,tF)
@@ -77,7 +74,7 @@ for (uh_tn, tn) in sol_t
   global _t_n
   _t_n += dt
   e = u(tn) - uh_tn
-  el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
+  el2 = sqrt(sum( ∫(l2(e))dΩ ))
   @test el2 < tol
 end
 
diff --git a/test/TransientFEsTests/runtests.jl b/test/TransientFEsTests/runtests.jl
index 05797ad..aec4392 100644
--- a/test/TransientFEsTests/runtests.jl
+++ b/test/TransientFEsTests/runtests.jl
@@ -26,4 +26,6 @@ using Test
 
 @testset "HeatEquationAutoDiffTests" begin include("HeatEquationAutoDiffTests.jl") end
 
+@testset "ForwardEulerHeatEquationTests" begin include("ForwardEulerHeatEquationTests.jl") end
+
 end # module