diff --git a/examples/03_solvers.jl b/examples/03_solvers.jl
index 7ece30a..b31929a 100644
--- a/examples/03_solvers.jl
+++ b/examples/03_solvers.jl
@@ -26,6 +26,10 @@ function ode_solver(
   return to_itensor(uₜ), nothing
 end
 
+function ode_solver(f⃗, H⃗₀, time_step, ψ₀; kwargs...)
+  return ode_solver(-im * TimeDependentSum(f⃗, H⃗₀), time_step, ψ₀; kwargs...)
+end
+
 function krylov_solver(
   H::TimeDependentSum, time_step, ψ₀; current_time=zero(time_step), outputlevel=0, kwargs...
 )
@@ -35,3 +39,7 @@ function krylov_solver(
   ψₜ, info = exponentiate(H(current_time), time_step, ψ₀; kwargs...)
   return ψₜ, info
 end
+
+function krylov_solver(f⃗, H⃗₀, time_step, ψ₀; kwargs...)
+  return krylov_solver(-im * TimeDependentSum(f⃗, H⃗₀), time_step, ψ₀; kwargs...)
+end
diff --git a/examples/03_tdvp_time_dependent.jl b/examples/03_tdvp_time_dependent.jl
index e012219..21d86a5 100644
--- a/examples/03_tdvp_time_dependent.jl
+++ b/examples/03_tdvp_time_dependent.jl
@@ -21,7 +21,7 @@ function main()
   # Set to 2 to get information about each bond/site
   # evolution, and 3 to get information about the
   # solver.
-  outputlevel = 1
+  outputlevel = 3
 
   # Frequency of time dependent terms
   ω₁ = 0.1
@@ -92,19 +92,12 @@ function main()
   println("#"^100)
   println()
 
-  function ode_solver(H⃗₀, time_step, ψ₀; kwargs...)
-    return ode_solver(
-      -im * TimeDependentSum(f⃗, H⃗₀),
-      time_step,
-      ψ₀;
-      solver_alg=ode_alg,
-      ode_kwargs...,
-      kwargs...,
-    )
+  function ode_solver_f⃗(H⃗₀, time_step, ψ₀; kwargs...)
+    return ode_solver(f⃗, H⃗₀, time_step, ψ₀; solver_alg=ode_alg, ode_kwargs..., kwargs...)
   end
 
   ψₜ_ode = tdvp(
-    ode_solver, H⃗₀, time_stop, ψ₀; time_step, maxdim, cutoff, nsite, outputlevel
+    ode_solver_f⃗, H⃗₀, time_stop, ψ₀; time_step, maxdim, cutoff, nsite, outputlevel
   )
 
   println()
@@ -117,13 +110,13 @@ function main()
   println("#"^100)
   println()
 
-  function krylov_solver(H⃗₀, time_step, ψ₀; kwargs...)
-    return krylov_solver(
-      -im * TimeDependentSum(f⃗, H⃗₀), time_step, ψ₀; krylov_kwargs..., kwargs...
-    )
+  function krylov_solver_f⃗(H⃗₀, time_step, ψ₀; kwargs...)
+    return krylov_solver(f⃗, H⃗₀, time_step, ψ₀; krylov_kwargs..., kwargs...)
   end
 
-  ψₜ_krylov = tdvp(krylov_solver, H⃗₀, time_stop, ψ₀; time_step, cutoff, nsite, outputlevel)
+  ψₜ_krylov = tdvp(
+    krylov_solver_f⃗, H⃗₀, time_stop, ψ₀; time_step, cutoff, nsite, outputlevel
+  )
 
   println()
   println("Finished running TDVP with Krylov solver")
@@ -136,7 +129,7 @@ function main()
   println()
 
   @disable_warn_order begin
-    ψₜ_full, _ = ode_solver(prod.(H⃗₀), time_stop, prod(ψ₀); outputlevel)
+    ψₜ_full, _ = ode_solver(f⃗, prod.(H⃗₀), time_stop, prod(ψ₀); outputlevel)
   end
 
   println()
diff --git a/src/tdvp.jl b/src/tdvp.jl
index 7133417..7c40795 100644
--- a/src/tdvp.jl
+++ b/src/tdvp.jl
@@ -1,8 +1,4 @@
-using ITensors:
-  Algorithm,
-  MPO,
-  MPS,
-  @Algorithm_str
+using ITensors: Algorithm, MPO, MPS, @Algorithm_str
 
 # Select solver function
 solver_function(solver_backend::String) = solver_function(Algorithm(solver_backend))
diff --git a/test/Project.toml b/test/Project.toml
index 0cdcad6..bfbc74b 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,6 +4,7 @@ ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/tdvp_time_dependent.jl b/test/tdvp_time_dependent.jl
index 4177673..eaa3d1b 100644
--- a/test/tdvp_time_dependent.jl
+++ b/test/tdvp_time_dependent.jl
@@ -1,43 +1,37 @@
 @eval module $(gensym())
 # These tests take too long to compile, skip for now.
-using ITensors: MPO, MPS, siteinds
+using ITensors: ITensors, MPO, MPS, siteinds
 using ITensorTDVP: ITensorTDVP, tdvp
+using LinearAlgebra: norm
 using Test: @test, @testset
 
 include(joinpath(pkgdir(ITensorTDVP), "examples", "03_models.jl"))
 include(joinpath(pkgdir(ITensorTDVP), "examples", "03_solvers.jl"))
 
-@testset "Time dependent Hamiltonian (eltype=$elt)" for elt in (
-  Float32, Float64, Complex{Float32}, Complex{Float64}
-)
-  function ode_solver(H⃗₀, time_step, ψ₀; kwargs...)
-    ω₁ = typeof(time_step)(0.1)
-    ω₂ = typeof(time_step)(0.2)
-    ω⃗ = [ω₁, ω₂]
+@testset "Time dependent Hamiltonian (eltype=$elt, conserve_qns=$conserve_qns)" for elt in
+                                                                                    (
+    Float32, Float64, Complex{Float32}, Complex{Float64}
+  ),
+  conserve_qns in [false, true]
+
+  # Frequencies
+  ω₁ = 0.1
+  ω₂ = 0.2
+
+  function ode_solver_f⃗(H⃗₀, time_step, ψ₀; kwargs...)
+    ω⃗ = typeof(time_step)[ω₁, ω₂]
     f⃗ = [t -> cos(ω * t) for ω in ω⃗]
-    ode_alg = Tsit5()
     tol = √eps(real(time_step))
-    ode_kwargs = (; reltol=tol, abstol=tol)
     return ode_solver(
-      -im * TimeDependentSum(f⃗, H⃗₀),
-      time_step,
-      ψ₀;
-      solver_alg=ode_alg,
-      ode_kwargs...,
-      kwargs...,
+      f⃗, H⃗₀, time_step, ψ₀; solver_alg=Tsit5(), reltol=tol, abstol=tol, kwargs...
     )
   end
 
-  function krylov_solver(H⃗₀, time_step, ψ₀; kwargs...)
-    ω₁ = typeof(time_step)(0.1)
-    ω₂ = typeof(time_step)(0.2)
-    ω⃗ = [ω₁, ω₂]
+  function krylov_solver_f⃗(H⃗₀, time_step, ψ₀; kwargs...)
+    ω⃗ = typeof(time_step)[ω₁, ω₂]
     f⃗ = [t -> cos(ω * t) for ω in ω⃗]
     tol = √eps(real(time_step))
-    krylov_kwargs = (; tol, eager=true)
-    return krylov_solver(
-      -im * TimeDependentSum(f⃗, H⃗₀), time_step, ψ₀; krylov_kwargs..., kwargs...
-    )
+    return krylov_solver(f⃗, H⃗₀, time_step, ψ₀; tol, eager=true, kwargs...)
   end
 
   n = 4
@@ -49,14 +43,14 @@ include(joinpath(pkgdir(ITensorTDVP), "examples", "03_solvers.jl"))
   maxdim = 100
   cutoff = √(eps(real(elt)))
   s = siteinds("S=1/2", n)
-  ℋ₁₀ = heisenberg(n; J=J₁, J2=elt(0))
-  ℋ₂₀ = heisenberg(n; J=elt(0), J2=J₂)
+  ℋ₁₀ = heisenberg(n; J=J₁, J2=zero(elt))
+  ℋ₂₀ = heisenberg(n; J=zero(elt), J2=J₂)
   ℋ⃗₀ = [ℋ₁₀, ℋ₂₀]
   H⃗₀ = [MPO(elt, ℋ₀, s) for ℋ₀ in ℋ⃗₀]
   ψ₀ = complex.(MPS(elt, s, j -> isodd(j) ? "↑" : "↓"))
-  ψₜ_ode = tdvp(ode_solver, H⃗₀, time_stop, ψ₀; time_step, maxdim, cutoff, nsite)
-  ψₜ_krylov = tdvp(krylov_solver, H⃗₀, time_stop, ψ₀; time_step, cutoff, nsite)
-  ψₜ_full, _ = ode_solver(prod.(H⃗₀), time_stop, prod(ψ₀))
+  ψₜ_ode = tdvp(ode_solver_f⃗, H⃗₀, time_stop, ψ₀; time_step, maxdim, cutoff, nsite)
+  ψₜ_krylov = tdvp(krylov_solver_f⃗, H⃗₀, time_stop, ψ₀; time_step, cutoff, nsite)
+  ψₜ_full, _ = ode_solver_f⃗(prod.(H⃗₀), time_stop, prod(ψ₀))
 
   @test ITensors.scalartype(ψ₀) == complex(elt)
   @test ITensors.scalartype(ψₜ_ode) == complex(elt)
diff --git a/test/test_examples.jl b/test/test_examples.jl
index 9592425..025b888 100644
--- a/test/test_examples.jl
+++ b/test/test_examples.jl
@@ -6,7 +6,7 @@ using Test: @testset
   examples_files = [
     "01_tdvp.jl",
     "02_dmrg-x.jl",
-    # "03_tdvp_time_dependent.jl",
+    "03_tdvp_time_dependent.jl",
     "04_tdvp_observers.jl",
     "05_tdvp_nonuniform_timesteps.jl",
   ]
diff --git a/test/test_tdvp_time_dependent.jl b/test/test_tdvp_time_dependent.jl
index 2478f8f..630f151 100644
--- a/test/test_tdvp_time_dependent.jl
+++ b/test/test_tdvp_time_dependent.jl
@@ -42,8 +42,7 @@ using Test: @test, @test_skip, @testset
     @test Hψ ≈ sum(i -> α * f⃗ₜ[i](t₀) * H⃗ᵣ[i](ψᵣ), eachindex(H⃗))
   end
   @testset "Time dependent TDVP" begin
-    # These tests take too long to compile, skip for now.
-    @test_skip include("tdvp_time_dependent.jl")
+    include("tdvp_time_dependent.jl")
   end
 end
 end