Not able to run the ODE-solvers properly

## Problem statement
Hi Gridap! Thanks for making a absolutely great FEM software library.

I am working on time-dependent problems, but unfortunately am I not able to get your solvers to work properly. Seems like `ThetaMethod` and `RungeKutta_BE_1_0_1` works as expected, but `RungeKutta_SDIRK_2_1_2` starts to diverge and `RungeKutta_TRBDF2_3_3_2` is crashing. For `Newmark` and `GeneralizedAlpha` am I not even able to compile. I am not sure if this is a user error or a bug.

Any idea about this?

## Minimalistic Example
I constructed a constructed a minimalistic example of the heat equation with Neumann condtions to test the solvers invidually.
```
using Dates
using Gridap

function run(;n, u_ex, dt, solver_choice="ThetaMethod", vtkdirname=nothing)
# Manufactured solution
f(t) = x -> ∂t(u_ex)(x,t) - Δ(u_ex(t))(x)
∇u_ex(t) = x -> ∇(u_ex(t))(x)

# Setup
L= 1.0
domain = (0,L,0,L)
partition = (n,n)
model = CartesianDiscreteModel(domain, partition)
order = 2
reffe = ReferenceFE(lagrangian, Float64, order)
V = FESpace( model, reffe, conformity=:H1)
U = TrialFESpace(V)

Ω = Triangulation(model)
Γ = BoundaryTriangulation(model)
degree = 2*order
dΩ = Measure(Ω,degree)
dΓ = Measure(Γ,degree)
n_Γ = get_normal_vector(Γ)

g(t) = ∇u_ex(t)⋅n_Γ
a(u,v) = ∫(∇(v)⋅∇(u))dΩ
b(v,t) = ∫(v*f(t))dΩ + ∫((g(t)⋅v))dΓ

res(t,u,v) = a(u,v) + ∫(∂t(u)*v)dΩ - b(v,t)
jac(t,u,du,v) = a(du,v)
jac_t(t,u,dut,v) = ∫(dut*v)dΩ

op = TransientFEOperator(res,jac,jac_t,U,V)

# Nonlinear/Linear solver
s =LUSolver()
# s = NLSolver(LUSolver();show_trace=true,method=:newton) #linesearch=BackTracking())

function choose_solver(solver_choice, solver_method, dt)
if solver_choice == "ThetaMethod"
return ThetaMethod(solver_method, dt, 1)
elseif solver_choice == "RungeKutta_BE_1_0_1"
return RungeKutta(solver_method, dt, :BE_1_0_1)
elseif solver_choice == "RungeKutta_SDIRK_2_1_2"
return RungeKutta(solver_method, dt, :SDIRK_2_1_2)
elseif solver_choice == "RungeKutta_TRBDF2_3_3_2"
return RungeKutta(solver_method, dt, :TRBDF2_3_3_2)
elseif solver_choice == "Newmark"
γ, β = 0.5, 0.25
return Newmark(solver_method,dt,γ,β)
elseif solver_choice == "GeneralizedAlpha"
ρ∞ = 1.0 # Equivalent to Newmark(0.5, 0.25)
return GeneralizedAlpha(solver_method, dt, ρ∞) # Does not compile
else
error("Invalid solver_choice: $solver_choice")
end
end

ode_solver = choose_solver(solver_choice, s, dt)

# Inital condition
t_0 = 0.0 # Edit: fixed floating error
T = 1.0
U_0 = interpolate_everywhere(0, U(0.0))

#################

U_h_t = solve(ode_solver, op, U_0, t_0, T)

ts = Float64[]
el2_ts = Float64[]
eh1_ts = Float64[]
eh_energy_ts = Float64[]

solname = vtkdirname*"/sol_dt_$dt"*"_n_$n"
mkpath(solname)
createpvd(solname*".pvd") do pvd
for (U_h, t) in U_h_t
e = u_ex(t) - U_h
pvd[t] = createvtk(Ω, solname*"_$t"*".vtu",cellfields=["u_h"=>U_h,"e"=>e])
el2_t = sqrt(sum( ∫(e*e)dΩ ))
eh1_t = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩ ))
eh_energy = sqrt(sum( ∫(∇(e)⋅∇(e) )*dΩ ))
push!( ts, t)
push!( el2_ts, el2_t )
push!( eh1_ts, eh1_t )
push!( eh_energy_ts, eh_energy)
end
end

el2s_L2 = sqrt(sum(dt* e_ti^2 for e_ti in el2_ts))
eh1s_L2 = sqrt(sum(dt* e_ti^2 for e_ti in eh1_ts))
ehs_energy_L2 = sqrt(sum(dt* e_ti^2 for e_ti in eh_energy_ts))

el2s_inf = maximum(abs.(el2_ts))
eh1s_inf = maximum(abs.(eh1_ts))
ehs_energy_inf = maximum(abs.(eh_energy_ts))

println("\nTesting solver choice $solver_choice")
println("L2 Norms:")
println("el2s_L2: $(round(el2s_L2, digits=4))")
println("eh1s_L2: $(round(eh1s_L2, digits=4))")
println("ehs_energy_L2: $(round(ehs_energy_L2, digits=4))")

println("Infinity Norms:")
println("el2s_inf: $(round(el2s_inf, digits=4))")
println("eh1s_inf: $(round(eh1s_inf, digits=4))")
println("ehs_energy_inf: $(round(ehs_energy_inf, digits=4))")
end

function main()
dirname= "figures/heat_eq/sim_"*string(Dates.now())
println(dirname)
mkpath(dirname)

u_ex(x,t) = sin(t)*cos(2π*x[1])*sin(2π*x[2])
u_ex(t::Real) = x -> u_ex(x,t)
run(n=2^5, dt=2^-4, u_ex=u_ex, vtkdirname=dirname, solver_choice="ThetaMethod")
run(n=2^5, dt=2^-4, u_ex=u_ex, vtkdirname=dirname, solver_choice="RungeKutta_BE_1_0_1")
run(n=2^5, dt=2^-4, u_ex=u_ex, vtkdirname=dirname, solver_choice="RungeKutta_SDIRK_2_1_2")
#run(n=2^5, dt=2^-4, u_ex=u_ex, vtkdirname=dirname, solver_choice="RungeKutta_TRBDF2_3_3_2")
#run(n=2^5, dt=2^-4, u_ex=u_ex, vtkdirname=dirname, solver_choice="Newmark")
#run(n=2^5, dt=2^-4, u_ex=u_ex, vtkdirname=dirname, solver_choice="GeneralizedAlpha")
end

main()

```
Where I get the following output

```
Testing solver choice ThetaMethod
L2 Norms:
el2s_L2: 0.0001
eh1s_L2: 0.0036
ehs_energy_L2: 0.0036
Infinity Norms:
el2s_inf: 0.0002
eh1s_inf: 0.0056
ehs_energy_inf: 0.0056

Testing solver choice RungeKutta_BE_1_0_1
L2 Norms:
el2s_L2: 0.0001
eh1s_L2: 0.0036
ehs_energy_L2: 0.0036
Infinity Norms:
el2s_inf: 0.0002
eh1s_inf: 0.0056
ehs_energy_inf: 0.0056

Testing solver choice RungeKutta_SDIRK_2_1_2
L2 Norms:
el2s_L2: 3.591747234711817e40
eh1s_L2: 1.2390496633141898e41
ehs_energy_L2: 1.1858488884248314e41
Infinity Norms:
el2s_inf: 1.4366961322166212e41
eh1s_inf: 4.95618880074577e41
ehs_energy_inf: 4.7433860956400755e41

```

As you can see `ThetaMethod` and `RungeKutta_BE_1_0_1` works properly and `RungeKutta_SDIRK_2_1_2` is diverging.

# Other Error Messages

If I run `RungeKutta_TRBDF2_3_3_2` I get the following crash report.

```
ERROR: LoadError: MethodError: no method matching update!(::Gridap.ODEs.ODETools.RungeKuttaNonlinearOperator, ::Int64, ::Vector{Vector{Float64}}, ::Int64)
Closest candidates are:
update!(::Gridap.ODEs.ODETools.RungeKuttaNonlinearOperator, ::Float64, ::AbstractVector, ::Int64) at ~/.julia/packages/Gridap/971dU/src/ODEs/ODETools/RungeKutta.jl:220
Stacktrace:
[1] solve_step!(uf::Vector{Float64}, solver::RungeKutta, op::Gridap.ODEs.TransientFETools.ODEOpFromFEOp{Gridap.ODEs.ODETools.Nonlinear}, u0::Vector{Float64}, t0::Int64, cache::Nothing)
@ Gridap.ODEs.ODETools ~/.julia/packages/Gridap/971dU/src/ODEs/ODETools/RungeKutta.jl:114
[2] solve_step!(uF::Vector{Float64}, solver::RungeKutta, op::Gridap.ODEs.TransientFETools.ODEOpFromFEOp{Gridap.ODEs.ODETools.Nonlinear}, u0::Vector{Float64}, t0::Int64)
@ Gridap.ODEs.ODETools ~/.julia/packages/Gridap/971dU/src/ODEs/ODETools/ODESolvers.jl:27
[3] iterate(sol::Gridap.ODEs.ODETools.GenericODESolution{Vector{Float64}})
@ Gridap.ODEs.ODETools ~/.julia/packages/Gridap/971dU/src/ODEs/ODETools/ODESolutions.jl:47
[4] iterate(sol::Gridap.ODEs.TransientFETools.TransientFESolution)
@ Gridap.ODEs.TransientFETools ~/.julia/packages/Gridap/971dU/src/ODEs/TransientFETools/TransientFESolutions.jl:65
[5] (::var"#131#146"{var"#u_ex#151", String, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Gridap.ODEs.TransientFETools.TransientFESolution, Gridap.CellData.GenericMeasure, Gridap.Geometry.BodyFittedTriangulation{2, 2, CartesianDiscreteModel{2, Float64, typeof(identity)}, CartesianGrid{2, Float64, typeof(identity)}, Gridap.Arrays.IdentityVector{Int64}}})(pvd::WriteVTK.CollectionFile)
@ Main ~/code/heat_eq_minimalistic.jl:80
[6] paraview_collection(f::var"#131#146"{var"#u_ex#151", String, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Gridap.ODEs.TransientFETools.TransientFESolution, Gridap.CellData.GenericMeasure, Gridap.Geometry.BodyFittedTriangulation{2, 2, CartesianDiscreteModel{2, Float64, typeof(identity)}, CartesianGrid{2, Float64, typeof(identity)}, Gridap.Arrays.IdentityVector{Int64}}}, args::String; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ WriteVTK ~/.julia/packages/WriteVTK/UoUwH/src/WriteVTK.jl:168
[7] paraview_collection(f::Function, args::String)
@ WriteVTK ~/.julia/packages/WriteVTK/UoUwH/src/WriteVTK.jl:164
[8] #createpvd#26
@ ~/.julia/packages/Gridap/971dU/src/Visualization/Vtk.jl:23 [inlined]
[9] createpvd
@ ~/.julia/packages/Gridap/971dU/src/Visualization/Vtk.jl:22 [inlined]
[10] run(; n::Int64, u_ex::var"#u_ex#151", dt::Float64, solver_choice::String, vtkdirname::String)
@ Main ~/code/heat_eq_minimalistic.jl:79
[11] main()
@ Main ~/code/heat_eq_minimalistic.jl:124
[12] top-level scope
@ ~/code/heat_eq_minimalistic.jl:130
[13] include(fname::String)
@ Base.MainInclude ./client.jl:476
[14] top-level scope
@ REPL[1]:1
in expression starting at /root/code/heat_eq_minimalistic.jl:130

```
And for both `Newman` and `GeneralizedAlpha` I simply get this error.
```
ERROR: LoadError: This function belongs to an interface definition and cannot be used.
```