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 8a1a629..6705884 100644
--- a/src/TransientFETools/TransientFEOperators.jl
+++ b/src/TransientFETools/TransientFEOperators.jl
@@ -193,11 +193,18 @@ function TransientFEOperator(res::Function,jac::Function,jac_t::Function,
 end
 
 function TransientFEOperator(res::Function,trial,test;order::Integer=1)
-  jacs = ()
-  for i in 1:order+1
+  function jac_0(t,x,dx0,dv)
+    function res_0(y)
+      x0 = (y,x[2:end]...)
+      res(t,x0,dv)
+    end
+    jacobian(res_0,x[1])
+  end
+  jacs = (jac_0,)
+  for i in 2:order+1
     function jac_i(t,x,dxi,dv)
       function res_i(y)
-        x[i] = y
+        x = (x[1:i-1]...,y,x[i+1:end]...)
         res(t,x,dv)
       end
       jacobian(res_i,x[i])
@@ -282,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 2de154d..759f014 100644
--- a/test/TransientFEsTests/ForwardEulerHeatEquationTests.jl
+++ b/test/TransientFEsTests/ForwardEulerHeatEquationTests.jl
@@ -1,4 +1,4 @@
-module HeatEquationTests
+module ForwardEulerHeatEquationTests
 
 using Gridap
 using ForwardDiff
@@ -30,25 +30,28 @@ model = CartesianDiscreteModel(domain,partition)
 
 order = 2
 
+reffe = ReferenceFE(lagrangian,Float64,order)
 V0 = FESpace(
-  reffe=lagrangian, order=order, valuetype=Float64,
-  conformity=:H1, model=model, dirichlet_tags="boundary")
+  model,
+  reffe,
+  conformity=:H1,
+  dirichlet_tags="boundary"
+)
 U = TransientTrialFESpace(V0,u)
 
-trian = Triangulation(model)
+Ω = Triangulation(model)
 degree = 2*order
-quad = CellQuadrature(trian,degree)
+dΩ = Measure(Ω,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
+res(t,(u,ut),v) = ∫( a(u,v) + ut*v - b(v,t) )dΩ
+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(U,V0,t_Ω)
+op = TransientFEOperator(res,jac,jac_t,U,V0)
 
 t0 = 0.0
 tF = 1.0
@@ -58,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)
@@ -73,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/HeatEquationAutoDiffTests.jl b/test/TransientFEsTests/HeatEquationAutoDiffTests.jl
index 9fb494c..06bdafe 100644
--- a/test/TransientFEsTests/HeatEquationAutoDiffTests.jl
+++ b/test/TransientFEsTests/HeatEquationAutoDiffTests.jl
@@ -50,11 +50,11 @@ dv = get_fe_basis(V0)
 du = get_trial_fe_basis(U₀)
 uh = FEFunction(U₀,rand(num_free_dofs(U₀)))
 
-cell_j = get_array(jac(0.5,uh,uh,du,dv))
-cell_j_t = get_array(jac_t(0.5,uh,uh,du,dv))
+cell_j = get_array(jac(0.5,(uh,uh),du,dv))
+cell_j_t = get_array(jac_t(0.5,(uh,uh),du,dv))
 
-cell_j_auto = get_array(jacobian(x->res(0.5,x,uh,dv),uh))
-cell_j_t_auto = get_array(jacobian(x->res(0.5,uh,x,dv),uh))
+cell_j_auto = get_array(jacobian(x->res(0.5,(x,uh),dv),uh))
+cell_j_t_auto = get_array(jacobian(x->res(0.5,(uh,x),dv),uh))
 
 test_array(cell_j_auto,cell_j,≈)
 test_array(cell_j_t_auto,cell_j_t,≈)
diff --git a/test/TransientFEsTests/StokesEquationAutoDiffTests.jl b/test/TransientFEsTests/StokesEquationAutoDiffTests.jl
index 487d6d6..1876df5 100644
--- a/test/TransientFEsTests/StokesEquationAutoDiffTests.jl
+++ b/test/TransientFEsTests/StokesEquationAutoDiffTests.jl
@@ -1,4 +1,4 @@
-module StokesEquationTests
+module StokesEquationAutoDiffTests
 
 using Gridap
 using ForwardDiff
@@ -6,7 +6,7 @@ using LinearAlgebra
 using Test
 using GridapODEs.ODETools
 using GridapODEs.TransientFETools
-using Gridap.FESpaces: get_algebraic_operator
+using Gridap.FESpaces
 using Gridap.Arrays: test_array
 
 # using GridapODEs.ODETools: ThetaMethodLinear
@@ -76,14 +76,18 @@ dy = get_fe_basis(Y)
 dx = get_trial_fe_basis(X₀)
 xh = FEFunction(X₀,rand(num_free_dofs(X₀)))
 
-cell_j = get_array(jac(0.5,xh,xh,dx,dy))
-cell_j_t = get_array(jac_t(0.5,xh,xh,dx,dy))
+cell_j = get_array(jac(0.5,(xh,xh),dx,dy))
+cell_j_t = get_array(jac_t(0.5,(xh,xh),dx,dy))
 
-cell_j_auto = get_array(jacobian(x->res(0.5,x,xh,dy),xh))
-cell_j_t_auto = get_array(jacobian(x->res(0.5,xh,x,dy),xh))
+cell_j_auto = get_array(jacobian(x->res(0.5,(x,xh),dy),xh))
+cell_j_t_auto = get_array(jacobian(x->res(0.5,(xh,x),dy),xh))
 
-test_array(cell_j_auto,cell_j,≈)
-test_array(cell_j_t_auto,cell_j_t,≈)
+for i in 1:length(cell_j)
+  test_array(cell_j[i].array[1,1],cell_j_auto[i].array[1,1],≈)
+  test_array(cell_j[i].array[1,2],cell_j_auto[i].array[1,2],≈)
+  test_array(cell_j[i].array[2,1],cell_j_auto[i].array[2,1],≈)
+  test_array(cell_j_t[i].array[1,1],cell_j_t_auto[i].array[1,1],≈)
+end
 
 op = TransientFEOperator(res,X,Y)
 
diff --git a/test/TransientFEsTests/runtests.jl b/test/TransientFEsTests/runtests.jl
index 95a0ca2..80dc0b1 100644
--- a/test/TransientFEsTests/runtests.jl
+++ b/test/TransientFEsTests/runtests.jl
@@ -24,4 +24,10 @@ using Test
 
 @testset "DGHeatEquationTests" begin include("DGHeatEquationTests.jl") end
 
+@testset "HeatEquationAutoDiffTests" begin include("HeatEquationAutoDiffTests.jl") end
+
+@testset "StokesEquationAutoDiffTests" begin include("StokesEquationAutoDiffTests.jl") end
+
+@testset "ForwardEulerHeatEquationTests" begin include("ForwardEulerHeatEquationTests.jl") end
+
 end # module