diff --git a/README.md b/README.md
index 3837353..74c4866 100644
--- a/README.md
+++ b/README.md
@@ -32,6 +32,6 @@ pkg> add ITensorTDVP
 
 ## About
 
-This package is effectively a generalization of the DMRG code in [ITensors.jl](https://github.com/ITensor/ITensors.jl), using the MPS/MPO types from that package. It provides a general MPS "solver" interface which allows us to implement other MPS/MPO optimization/solver functionality like DMRG (`ITensorTDVP.dmrg` or `eigsolve`), TDVP (`tdvp` or `exponentiate`), linear solving (`linsolve`), DMRG-X (`dmrg_x`), etc. while sharing most of the code across those different functions. Therefore, it effectively supercedes the DMRG functionality in ITensors.jl (`dmrg`), and provides its own `ITensorTDVP.dmrg`/`eigsolve` function that is essentially the same as the `dmrg` function from ITensors.jl. This package is fairly stable and appropriate for general use. The primary missing feature is a lack of modern subspace expansion tools for methods like TDVP. However, 2-site TDVP or TEBD is often sufficient for performing subspace expansion (except when [it's not](https://arxiv.org/abs/2005.06104)).
+This package is effectively a generalization of the DMRG code in [ITensors.jl](https://github.com/ITensor/ITensors.jl), using the MPS/MPO types from that package. It provides a general MPS "solver" interface which allows us to implement other MPS/MPO optimization/solver functionality like DMRG (`ITensorTDVP.dmrg`), TDVP (`ITensorTDVP.tdvp`), linear solving (`ITensorTDVP.linsolve`/`KrylovKit.linsolve`), DMRG-X (`ITensorTDVP.dmrg_x`), etc. while sharing most of the code across those different functions. Therefore, it effectively supercedes the DMRG functionality in ITensors.jl (`dmrg`), and provides its own `ITensorTDVP.dmrg` function that is essentially the same as the `dmrg` function from ITensors.jl. This package is fairly stable and appropriate for general use. The primary missing feature is a lack of modern subspace expansion tools for methods like TDVP. However, 2-site TDVP or TEBD is often sufficient for performing subspace expansion (except when [it's not](https://arxiv.org/abs/2005.06104)).
 
 However, note that future developments, including modern subspace expansion tools, are being developed in our next-generation tensor network library [ITensorNetworks.jl](https://github.com/mtfishman/ITensorNetworks.jl). The goal of that package is to provide contraction, optimization, and evolution tools for general tensor networks, as well as methods like DMRG, TDVP, and linear solving for tree tensor networks, and the eventual goal is to replace this package which is limited to solvers for just MPS/MPO (linear/path graph) tensor networks. However, ITensorNetworks.jl is under heavy development and is _not_ meant for general usage at the moment, except for those who are brave enough to handle missing features and breaking interfaces. Basically, for the average user who wants stable and reliable code, if you need to use MPS-based TDVP or linear solving, you should use this package for the time being.
diff --git a/examples/01_tdvp.jl b/examples/01_tdvp.jl
index 0fad273..c386c0a 100644
--- a/examples/01_tdvp.jl
+++ b/examples/01_tdvp.jl
@@ -1,42 +1,47 @@
-using ITensors
-using ITensorTDVP
-
-n = 10
-s = siteinds("S=1/2", n)
-
-function heisenberg(n)
-  os = OpSum()
-  for j in 1:(n - 1)
-    os += 0.5, "S+", j, "S-", j + 1
-    os += 0.5, "S-", j, "S+", j + 1
-    os += "Sz", j, "Sz", j + 1
+using ITensors: ITensors, MPO, OpSum, inner, randomMPS, siteinds
+using ITensorTDVP: ITensorTDVP, tdvp
+
+function main()
+  n = 10
+  s = siteinds("S=1/2", n)
+
+  function heisenberg(n)
+    os = OpSum()
+    for j in 1:(n - 1)
+      os += 0.5, "S+", j, "S-", j + 1
+      os += 0.5, "S-", j, "S+", j + 1
+      os += "Sz", j, "Sz", j + 1
+    end
+    return os
   end
-  return os
-end
 
-H = MPO(heisenberg(n), s)
-ψ = randomMPS(s, "↑"; linkdims=10)
+  H = MPO(heisenberg(n), s)
+  ψ = randomMPS(s, "↑"; linkdims=10)
+
+  @show inner(ψ', H, ψ) / inner(ψ, ψ)
 
-@show inner(ψ', H, ψ) / inner(ψ, ψ)
+  ϕ = tdvp(
+    H,
+    -1.0,
+    ψ;
+    nsweeps=20,
+    reverse_step=false,
+    normalize=true,
+    maxdim=30,
+    cutoff=1e-10,
+    outputlevel=1,
+  )
 
-ϕ = tdvp(
-  H,
-  -1.0,
-  ψ;
-  nsweeps=20,
-  reverse_step=false,
-  normalize=true,
-  maxdim=30,
-  cutoff=1e-10,
-  outputlevel=1,
-)
+  @show inner(ϕ', H, ϕ) / inner(ϕ, ϕ)
 
-@show inner(ϕ', H, ϕ) / inner(ϕ, ϕ)
+  e2, ϕ2 = ITensors.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10)
 
-e2, ϕ2 = dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10)
+  @show inner(ϕ2', H, ϕ2) / inner(ϕ2, ϕ2)
 
-@show inner(ϕ2', H, ϕ2) / inner(ϕ2, ϕ2)
+  ϕ3 = ITensorTDVP.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10, outputlevel=1)
 
-ϕ3 = ITensorTDVP.dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10, outputlevel=1)
+  @show inner(ϕ3', H, ϕ3) / inner(ϕ3, ϕ3)
+  return nothing
+end
 
-@show inner(ϕ3', H, ϕ3) / inner(ϕ3, ϕ3)
+main()
diff --git a/examples/02_dmrg-x.jl b/examples/02_dmrg-x.jl
index 2d6ce60..6d9e670 100644
--- a/examples/02_dmrg-x.jl
+++ b/examples/02_dmrg-x.jl
@@ -1,44 +1,48 @@
-using ITensors
-using ITensorTDVP
-using LinearAlgebra
-
-function heisenberg(n; h=zeros(n))
-  os = OpSum()
-  for j in 1:(n - 1)
-    os += 0.5, "S+", j, "S-", j + 1
-    os += 0.5, "S-", j, "S+", j + 1
-    os += "Sz", j, "Sz", j + 1
-  end
-  for j in 1:n
-    if h[j] ≠ 0
-      os -= h[j], "Sz", j
+using ITensors: MPO, MPS, OpSum, inner, siteinds
+using ITensorTDVP: dmrg_x
+using Random: Random
+
+function main()
+  function heisenberg(n; h=zeros(n))
+    os = OpSum()
+    for j in 1:(n - 1)
+      os += 0.5, "S+", j, "S-", j + 1
+      os += 0.5, "S-", j, "S+", j + 1
+      os += "Sz", j, "Sz", j + 1
+    end
+    for j in 1:n
+      if h[j] ≠ 0
+        os -= h[j], "Sz", j
+      end
     end
+    return os
   end
-  return os
-end
 
-n = 10
-s = siteinds("S=1/2", n)
+  n = 10
+  s = siteinds("S=1/2", n)
 
-using Random
-Random.seed!(12)
+  Random.seed!(12)
 
-# MBL when W > 3.5-4
-W = 12
-# Random fields h ∈ [-W, W]
-h = W * (2 * rand(n) .- 1)
-H = MPO(heisenberg(n; h), s)
+  # MBL when W > 3.5-4
+  W = 12
+  # Random fields h ∈ [-W, W]
+  h = W * (2 * rand(n) .- 1)
+  H = MPO(heisenberg(n; h), s)
 
-initstate = rand(["↑", "↓"], n)
-ψ = MPS(s, initstate)
+  initstate = rand(["↑", "↓"], n)
+  ψ = MPS(s, initstate)
 
-dmrg_x_kwargs = (
-  nsweeps=10, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=1
-)
+  dmrg_x_kwargs = (
+    nsweeps=10, reverse_step=false, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=1
+  )
 
-ϕ = dmrg_x(H, ψ; dmrg_x_kwargs...)
+  ϕ = dmrg_x(H, ψ; dmrg_x_kwargs...)
+
+  @show inner(ψ', H, ψ) / inner(ψ, ψ)
+  @show inner(H, ψ, H, ψ) - inner(ψ', H, ψ)^2
+  @show inner(ϕ', H, ϕ) / inner(ϕ, ϕ)
+  @show inner(H, ϕ, H, ϕ) - inner(ϕ', H, ϕ)^2
+  return nothing
+end
 
-@show inner(ψ', H, ψ) / inner(ψ, ψ)
-@show inner(H, ψ, H, ψ) - inner(ψ', H, ψ)^2
-@show inner(ϕ', H, ϕ) / inner(ϕ, ϕ)
-@show inner(H, ϕ, H, ϕ) - inner(ϕ', H, ϕ)^2
+main()
diff --git a/examples/03_tdvp_time_dependent.jl b/examples/03_tdvp_time_dependent.jl
index 124caf3..e012219 100644
--- a/examples/03_tdvp_time_dependent.jl
+++ b/examples/03_tdvp_time_dependent.jl
@@ -149,6 +149,7 @@ function main()
 
   @show 1 - abs(inner(prod(ψₜ_ode), ψₜ_full))
   @show 1 - abs(inner(prod(ψₜ_krylov), ψₜ_full))
+  return nothing
 end
 
 main()
diff --git a/examples/04_tdvp_observers.jl b/examples/04_tdvp_observers.jl
index 8ab3f8a..f2c056d 100644
--- a/examples/04_tdvp_observers.jl
+++ b/examples/04_tdvp_observers.jl
@@ -1,81 +1,86 @@
-using ITensors
-using ITensorTDVP
-using Observers
+using ITensors: MPO, MPS, OpSum, expect, siteinds
+using ITensorTDVP: tdvp
+using Observers: observer
 
-function heisenberg(N)
-  os = OpSum()
-  for j in 1:(N - 1)
-    os += 0.5, "S+", j, "S-", j + 1
-    os += 0.5, "S-", j, "S+", j + 1
-    os += "Sz", j, "Sz", j + 1
+function main()
+  function heisenberg(N)
+    os = OpSum()
+    for j in 1:(N - 1)
+      os += 0.5, "S+", j, "S-", j + 1
+      os += 0.5, "S-", j, "S+", j + 1
+      os += "Sz", j, "Sz", j + 1
+    end
+    return os
   end
-  return os
-end
 
-N = 10
-cutoff = 1e-12
-tau = 0.1
-ttotal = 1.0
+  N = 10
+  cutoff = 1e-12
+  tau = 0.1
+  ttotal = 1.0
 
-s = siteinds("S=1/2", N; conserve_qns=true)
-H = MPO(heisenberg(N), s)
+  s = siteinds("S=1/2", N; conserve_qns=true)
+  H = MPO(heisenberg(N), s)
 
-function step(; sweep, bond, half_sweep)
-  if bond == 1 && half_sweep == 2
-    return sweep
+  function step(; sweep, bond, half_sweep)
+    if bond == 1 && half_sweep == 2
+      return sweep
+    end
+    return nothing
   end
-  return nothing
-end
 
-function current_time(; current_time, bond, half_sweep)
-  if bond == 1 && half_sweep == 2
-    return current_time
+  function current_time(; current_time, bond, half_sweep)
+    if bond == 1 && half_sweep == 2
+      return current_time
+    end
+    return nothing
   end
-  return nothing
-end
 
-function measure_sz(; psi, bond, half_sweep)
-  if bond == 1 && half_sweep == 2
-    return expect(psi, "Sz"; sites=N ÷ 2)
+  function measure_sz(; psi, bond, half_sweep)
+    if bond == 1 && half_sweep == 2
+      return expect(psi, "Sz"; sites=N ÷ 2)
+    end
+    return nothing
   end
-  return nothing
-end
 
-function return_state(; psi, bond, half_sweep)
-  if bond == 1 && half_sweep == 2
-    return psi
+  function return_state(; psi, bond, half_sweep)
+    if bond == 1 && half_sweep == 2
+      return psi
+    end
+    return nothing
   end
-  return nothing
-end
 
-obs = observer(
-  "steps" => step, "times" => current_time, "psis" => return_state, "Sz" => measure_sz
-)
+  obs = observer(
+    "steps" => step, "times" => current_time, "psis" => return_state, "Sz" => measure_sz
+  )
 
-psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")
-psi_f = tdvp(
-  H,
-  -im * ttotal,
-  psi;
-  time_step=-im * tau,
-  cutoff,
-  outputlevel=1,
-  normalize=false,
-  (observer!)=obs,
-)
+  psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")
+  psi_f = tdvp(
+    H,
+    -im * ttotal,
+    psi;
+    time_step=-im * tau,
+    cutoff,
+    outputlevel=1,
+    normalize=false,
+    (observer!)=obs,
+  )
 
-steps = obs.steps
-times = obs.times
-psis = obs.psis
-Sz = obs.Sz
+  steps = obs.steps
+  times = obs.times
+  psis = obs.psis
+  Sz = obs.Sz
 
-println("\nResults")
-println("=======")
-for n in 1:length(steps)
-  print("step = ", steps[n])
-  print(", time = ", round(times[n]; digits=3))
-  print(", |⟨ψⁿ|ψⁱ⟩| = ", round(abs(inner(psis[n], psi)); digits=3))
-  print(", |⟨ψⁿ|ψᶠ⟩| = ", round(abs(inner(psis[n], psi_f)); digits=3))
-  print(", ⟨Sᶻ⟩ = ", round(Sz[n]; digits=3))
-  println()
+  println("\nResults")
+  println("=======")
+  for n in 1:length(steps)
+    print("step = ", steps[n])
+    print(", time = ", round(times[n]; digits=3))
+    print(", |⟨ψⁿ|ψⁱ⟩| = ", round(abs(inner(psis[n], psi)); digits=3))
+    print(", |⟨ψⁿ|ψᶠ⟩| = ", round(abs(inner(psis[n], psi_f)); digits=3))
+    print(", ⟨Sᶻ⟩ = ", round(Sz[n]; digits=3))
+    println()
+  end
+  return nothing
 end
+
+main()
diff --git a/examples/05_tdvp_nonuniform_timesteps.jl b/examples/05_tdvp_nonuniform_timesteps.jl
index bf6f67a..0a76206 100644
--- a/examples/05_tdvp_nonuniform_timesteps.jl
+++ b/examples/05_tdvp_nonuniform_timesteps.jl
@@ -1,46 +1,52 @@
-using ITensors
-using ITensorTDVP
+using ITensors: MPO, MPS, OpSum, ProjMPO, siteinds
+using KrylovKit: exponentiate
+using Observers: observer
 
 include("05_utils.jl")
 
-function heisenberg(N)
-  os = OpSum()
-  for j in 1:(N - 1)
-    os += 0.5, "S+", j, "S-", j + 1
-    os += 0.5, "S-", j, "S+", j + 1
-    os += "Sz", j, "Sz", j + 1
+function main()
+  function heisenberg(N)
+    os = OpSum()
+    for j in 1:(N - 1)
+      os += 0.5, "S+", j, "S-", j + 1
+      os += 0.5, "S-", j, "S+", j + 1
+      os += "Sz", j, "Sz", j + 1
+    end
+    return os
   end
-  return os
-end
 
-N = 10
-cutoff = 1e-12
-outputlevel = 1
-nsteps = 10
-time_steps = [n ≤ 2 ? -0.2im : -0.1im for n in 1:nsteps]
-
-obs = observer("times" => (; current_time) -> current_time, "psis" => (; psi) -> psi)
-
-s = siteinds("S=1/2", N; conserve_qns=true)
-H = MPO(heisenberg(N), s)
-
-psi0 = MPS(s, n -> isodd(n) ? "Up" : "Dn")
-psi = tdvp_nonuniform_timesteps(
-  ProjMPO(H), psi0; time_steps, cutoff, outputlevel, (step_observer!)=obs
-)
-
-times = obs.times
-psis = obs.psis
-
-println("\nResults")
-println("=======")
-print("step = ", 0)
-print(", time = ", zero(ComplexF64))
-print(", ⟨Sᶻ⟩ = ", round(expect(psi0, "Sz"; sites=N ÷ 2); digits=3))
-println()
-for n in 1:length(times)
-  print("step = ", n)
-  print(", time = ", round(times[n]; digits=3))
-  print(", ⟨Sᶻ⟩ = ", round(expect(psis[n], "Sz"; sites=N ÷ 2); digits=3))
+  N = 10
+  cutoff = 1e-12
+  outputlevel = 1
+  nsteps = 10
+  time_steps = [n ≤ 2 ? -0.2im : -0.1im for n in 1:nsteps]
+
+  obs = observer("times" => (; current_time) -> current_time, "psis" => (; psi) -> psi)
+
+  s = siteinds("S=1/2", N; conserve_qns=true)
+  H = MPO(heisenberg(N), s)
+
+  psi0 = MPS(s, n -> isodd(n) ? "Up" : "Dn")
+  psi = tdvp_nonuniform_timesteps(
+    ProjMPO(H), psi0; time_steps, cutoff, outputlevel, (step_observer!)=obs
+  )
+
+  times = obs.times
+  psis = obs.psis
+
+  println("\nResults")
+  println("=======")
+  print("step = ", 0)
+  print(", time = ", zero(ComplexF64))
+  print(", ⟨Sᶻ⟩ = ", round(expect(psi0, "Sz"; sites=N ÷ 2); digits=3))
   println()
+  for n in 1:length(times)
+    print("step = ", n)
+    print(", time = ", round(times[n]; digits=3))
+    print(", ⟨Sᶻ⟩ = ", round(expect(psis[n], "Sz"; sites=N ÷ 2); digits=3))
+    println()
+  end
+  return nothing
 end
+
+main()
diff --git a/examples/05_utils.jl b/examples/05_utils.jl
index b168b14..5d61efd 100644
--- a/examples/05_utils.jl
+++ b/examples/05_utils.jl
@@ -1,9 +1,7 @@
-using ITensors
-using ITensorTDVP
-using Observers
-using Printf
-
-using ITensorTDVP: tdvp_solver, tdvp_step, process_sweeps, TDVPOrder
+using ITensors: MPS, maxlinkdim
+using ITensorTDVP: TDVPOrder, process_sweeps, tdvp_solver, tdvp_step, process_sweeps
+using Observers: observer, update!
+using Printf: @printf
 
 function tdvp_nonuniform_timesteps(
   solver,
@@ -14,10 +12,15 @@ function tdvp_nonuniform_timesteps(
   time_start=0.0,
   order=2,
   (step_observer!)=observer(),
+  maxdim=ITensorTDVP.default_maxdim(),
+  mindim=ITensorTDVP.default_mindim(),
+  cutoff=ITensorTDVP.default_cutoff(),
+  noise=ITensorTDVP.default_noise(),
+  outputlevel=ITensorTDVP.default_outputlevel(),
   kwargs...,
 )
   nsweeps = length(time_steps)
-  maxdim, mindim, cutoff, noise = process_sweeps(; nsweeps, kwargs...)
+  maxdim, mindim, cutoff, noise = process_sweeps(; nsweeps, maxdim, mindim, cutoff, noise)
   tdvp_order = TDVPOrder(order, Base.Forward)
   current_time = time_start
   for sw in 1:nsweeps
@@ -55,6 +58,29 @@ function tdvp_nonuniform_timesteps(
   return psi
 end
 
-function tdvp_nonuniform_timesteps(H, psi::MPS; kwargs...)
-  return tdvp_nonuniform_timesteps(tdvp_solver(; kwargs...), H, psi; kwargs...)
+function tdvp_nonuniform_timesteps(
+  H,
+  psi::MPS;
+  ishermitian=ITensorTDVP.default_ishermitian(),
+  issymmetric=ITensorTDVP.default_issymmetric(),
+  solver_tol=ITensorTDVP.default_solver_tol(exponentiate),
+  solver_krylovdim=ITensorTDVP.default_solver_krylovdim(exponentiate),
+  solver_maxiter=ITensorTDVP.default_solver_maxiter(exponentiate),
+  solver_outputlevel=ITensorTDVP.default_solver_outputlevel(exponentiate),
+  kwargs...,
+)
+  return tdvp_nonuniform_timesteps(
+    tdvp_solver(
+      exponentiate;
+      ishermitian,
+      issymmetric,
+      solver_tol,
+      solver_krylovdim,
+      solver_maxiter,
+      solver_outputlevel,
+    ),
+    H,
+    psi;
+    kwargs...,
+  )
 end
diff --git a/src/tdvp_generic.jl b/src/tdvp_generic.jl
index 066e78f..6279ac8 100644
--- a/src/tdvp_generic.jl
+++ b/src/tdvp_generic.jl
@@ -7,6 +7,7 @@ using ITensors:
   checkdone!,
   disk,
   linkind,
+  maxlinkdim,
   permute
 
 function _tdvp_compute_nsweeps(t; time_step=default_time_step(t), nsweeps=default_nsweeps())
diff --git a/test/Project.toml b/test/Project.toml
index 7f0787b..0cdcad6 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -5,4 +5,5 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/test_examples.jl b/test/test_examples.jl
new file mode 100644
index 0000000..9592425
--- /dev/null
+++ b/test/test_examples.jl
@@ -0,0 +1,19 @@
+@eval module $(gensym())
+using Suppressor: @suppress
+using ITensorTDVP: ITensorTDVP
+using Test: @testset
+@testset "Run examples" begin
+  examples_files = [
+    "01_tdvp.jl",
+    "02_dmrg-x.jl",
+    # "03_tdvp_time_dependent.jl",
+    "04_tdvp_observers.jl",
+    "05_tdvp_nonuniform_timesteps.jl",
+  ]
+  examples_path = joinpath(pkgdir(ITensorTDVP), "examples")
+  @testset "Running example file $f" for f in examples_files
+    println("Running example file $f")
+    @suppress include(joinpath(examples_path, f))
+  end
+end
+end
diff --git a/test/test_tdvp.jl b/test/test_tdvp.jl
index 38bcb32..97c7146 100644
--- a/test/test_tdvp.jl
+++ b/test/test_tdvp.jl
@@ -1,6 +1,20 @@
 @eval module $(gensym())
 using ITensors:
-  ITensors, AbstractObserver, ITensor, MPO, MPS, OpSum, randomMPS, inner, op, siteinds
+  ITensors,
+  AbstractObserver,
+  ITensor,
+  MPO,
+  MPS,
+  OpSum,
+  apply,
+  dag,
+  expect,
+  inner,
+  noprime,
+  op,
+  prime,
+  randomMPS,
+  siteinds
 using ITensorTDVP: tdvp
 using KrylovKit: exponentiate
 using LinearAlgebra: norm