diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 333e3af4e1..1e363e5b08 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -192,7 +192,7 @@ const D = Differential(t)
 PrecompileTools.@compile_workload begin
     using ModelingToolkit
     @variables x(ModelingToolkit.t_nounits)
-    @named sys = ODESystem([ModelingToolkit.D_nounits(x) ~ -x])
+    @named sys = ODESystem([ModelingToolkit.D_nounits(x) ~ -x], ModelingToolkit.t_nounits)
     prob = ODEProblem(structural_simplify(sys), [x => 30.0], (0, 100), [], jac = true)
 end
 
diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index 880fdd1c6b..751d8f7ef7 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -230,7 +230,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
         metadata, gui_metadata, checks = checks)
 end
 
-function ODESystem(eqs, iv = nothing; kwargs...)
+function ODESystem(eqs, iv; kwargs...)
     eqs = scalarize(eqs)
     # NOTE: this assumes that the order of algebraic equations doesn't matter
     diffvars = OrderedSet()
diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index 3edf34252c..4267ca9d6d 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -1,4 +1,4 @@
-function System(eqs::AbstractVector{<:Equation}, iv = nothing, args...; name = nothing,
+function System(eqs::AbstractVector{<:Equation}, iv, args...; name = nothing,
         kw...)
     ODESystem(eqs, iv, args...; name, kw..., checks = false)
 end
diff --git a/test/clock.jl b/test/clock.jl
index c694340ec1..63904bd01d 100644
--- a/test/clock.jl
+++ b/test/clock.jl
@@ -22,7 +22,7 @@ eqs = [yd ~ Sample(t, dt)(y)
     u ~ Hold(ud)
     D(x) ~ -x + u
     y ~ x]
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 # compute equation and variables' time domains
 #TODO: test linearize
 
@@ -110,7 +110,7 @@ z′(k + 1) ~ z(k) + yd
 z(k + 1)  ~ z′(k)
 =#
 ]
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 ss = structural_simplify(sys);
 
 Tf = 1.0
@@ -174,7 +174,7 @@ eqs = [
     u ~ Hold(ud1) + Hold(ud2)
     D(x) ~ -x + u
     y ~ x]
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 ci, varmap = infer_clocks(sys)
 
 d = Clock(t, dt)
diff --git a/test/constants.jl b/test/constants.jl
index 373c37464a..2427638703 100644
--- a/test/constants.jl
+++ b/test/constants.jl
@@ -9,7 +9,7 @@ UMT = ModelingToolkit.UnitfulUnitCheck
 @variables t x(t) w(t)
 D = Differential(t)
 eqs = [D(x) ~ a]
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 prob = ODEProblem(complete(sys), [0], [0.0, 1.0], [])
 sol = solve(prob, Tsit5())
 
@@ -19,7 +19,7 @@ newsys = MT.eliminate_constants(sys)
 # Test structural_simplify substitutions & observed values
 eqs = [D(x) ~ 1,
     w ~ a]
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 # Now eliminate the constants first
 simp = structural_simplify(sys)
 @test equations(simp) == [D(x) ~ 1.0]
@@ -31,7 +31,7 @@ UMT.get_unit(β)
 @variables t [unit = u"s"] x(t) [unit = u"m"]
 D = Differential(t)
 eqs = [D(x) ~ β]
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 simp = structural_simplify(sys)
 
 @test isempty(MT.collect_constants(nothing))
diff --git a/test/dae_jacobian.jl b/test/dae_jacobian.jl
index de5a0542a5..8c68df767a 100644
--- a/test/dae_jacobian.jl
+++ b/test/dae_jacobian.jl
@@ -36,7 +36,7 @@ sol1 = solve(prob1, IDA(linear_solver = :KLU))
 eqs = [D(u1) ~ p1 * u1 - u1 * u2,
     D(u2) ~ u1 * u2 - p2 * u2]
 
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 
 u0 = [u1 => 1.0,
     u2 => 1.0]
diff --git a/test/distributed.jl b/test/distributed.jl
index 93f6aaebdf..3a53b9951e 100644
--- a/test/distributed.jl
+++ b/test/distributed.jl
@@ -13,7 +13,7 @@ addprocs(2)
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z]
 
-@everywhere @named de = ODESystem(eqs)
+@everywhere @named de = ODESystem(eqs, t)
 @everywhere de = complete(de)
 @everywhere ode_func = ODEFunction(de, [x, y, z], [σ, ρ, β])
 
diff --git a/test/dq_units.jl b/test/dq_units.jl
index 758b468e8b..c4069be17d 100644
--- a/test/dq_units.jl
+++ b/test/dq_units.jl
@@ -20,7 +20,7 @@ using ModelingToolkit: t, D
 eqs = [D(E) ~ P - E / τ
     0 ~ P]
 @test MT.validate(eqs)
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 
 @test !MT.validate(D(D(E)) ~ P)
 @test !MT.validate(0 ~ P + E * τ)
@@ -30,13 +30,13 @@ eqs = [D(E) ~ P - E / τ
 ODESystem(eqs, t, [E, P, t], [τ], name = :sys, checks = MT.CheckUnits)
 eqs = [D(E) ~ P - E / τ
     0 ~ P + E * τ]
-@test_throws MT.ValidationError ODESystem(eqs, name = :sys, checks = MT.CheckAll)
-@test_throws MT.ValidationError ODESystem(eqs, name = :sys, checks = true)
-ODESystem(eqs, name = :sys, checks = MT.CheckNone)
-ODESystem(eqs, name = :sys, checks = false)
-@test_throws MT.ValidationError ODESystem(eqs, name = :sys,
+@test_throws MT.ValidationError ODESystem(eqs, t, name = :sys, checks = MT.CheckAll)
+@test_throws MT.ValidationError ODESystem(eqs, t, name = :sys, checks = true)
+ODESystem(eqs, t, name = :sys, checks = MT.CheckNone)
+ODESystem(eqs, t, name = :sys, checks = false)
+@test_throws MT.ValidationError ODESystem(eqs, t, name = :sys,
     checks = MT.CheckComponents | MT.CheckUnits)
-@named sys = ODESystem(eqs, checks = MT.CheckComponents)
+@named sys = ODESystem(eqs, t, checks = MT.CheckComponents)
 @test_throws MT.ValidationError ODESystem(eqs, t, [E, P, t], [τ], name = :sys,
     checks = MT.CheckUnits)
 
@@ -69,7 +69,7 @@ good_eqs = [connect(p1, p2)]
 @variables x(t)[1:3] [unit = u"m"]
 @parameters v[1:3]=[1, 2, 3] [unit = u"m/s"]
 eqs = D.(x) .~ v
-ODESystem(eqs, name = :sys)
+ODESystem(eqs, t, name = :sys)
 
 # Nonlinear system
 @parameters a [unit = u"kg"^-1]
@@ -104,12 +104,12 @@ noiseeqs = [0.1u"W" 0.1u"W"
 @parameters v [unit = u"m/s"] r [unit = u"m"^3 / u"s"]
 eqs = [D(L) ~ v,
     V ~ L^3]
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 sys_simple = structural_simplify(sys)
 
 eqs = [D(V) ~ r,
     V ~ L^3]
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 sys_simple = structural_simplify(sys)
 
 @variables V [unit = u"m"^3] L [unit = u"m"]
diff --git a/test/extensions/bifurcationkit.jl b/test/extensions/bifurcationkit.jl
index a0e56f5028..ec5ec94604 100644
--- a/test/extensions/bifurcationkit.jl
+++ b/test/extensions/bifurcationkit.jl
@@ -60,7 +60,7 @@ let
     @variables x(t) y(t) z(t)
     eqs = [D(x) ~ -x + a * y + x^2 * y,
         D(y) ~ b - a * y - x^2 * y]
-    @named sys = ODESystem(eqs)
+    @named sys = ODESystem(eqs, t)
     sys = complete(sys)
     # Creates BifurcationProblem
     bprob = BifurcationProblem(sys,
diff --git a/test/input_output_handling.jl b/test/input_output_handling.jl
index f9f7db5453..89b072916f 100644
--- a/test/input_output_handling.jl
+++ b/test/input_output_handling.jl
@@ -158,7 +158,7 @@ eqs = [
     D(x) ~ -x + u,
 ]
 
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 f, dvs, ps = ModelingToolkit.generate_control_function(sys, simplify = true)
 
 @test isequal(dvs[], x)
diff --git a/test/labelledarrays.jl b/test/labelledarrays.jl
index 3195aa15b0..7808cf6f2e 100644
--- a/test/labelledarrays.jl
+++ b/test/labelledarrays.jl
@@ -12,7 +12,7 @@ eqs = [D(x) ~ σ * (y - x),
     D(y) ~ t * x * (ρ - z) - y,
     D(z) ~ x * y - β * z]
 
-@named de = ODESystem(eqs)
+@named de = ODESystem(eqs, t)
 de = complete(de)
 ff = ODEFunction(de, [x, y, z], [σ, ρ, β], jac = true)
 
diff --git a/test/linearity.jl b/test/linearity.jl
index 1c83b550b1..0293110c7f 100644
--- a/test/linearity.jl
+++ b/test/linearity.jl
@@ -11,16 +11,16 @@ eqs = [D(x) ~ σ * (y - x),
     D(y) ~ -z - y,
     D(z) ~ y - β * z]
 
-@test ModelingToolkit.islinear(@named sys = ODESystem(eqs))
+@test ModelingToolkit.islinear(@named sys = ODESystem(eqs, t))
 
 eqs2 = [D(x) ~ σ * (y - x),
     D(y) ~ -z - 1 / y,
     D(z) ~ y - β * z]
 
-@test !ModelingToolkit.islinear(@named sys = ODESystem(eqs2))
+@test !ModelingToolkit.islinear(@named sys = ODESystem(eqs2, t))
 
 eqs3 = [D(x) ~ σ * (y - x),
     D(y) ~ -z - y,
     D(z) ~ y - β * z + 1]
 
-@test ModelingToolkit.isaffine(@named sys = ODESystem(eqs))
+@test ModelingToolkit.isaffine(@named sys = ODESystem(eqs, t))
diff --git a/test/lowering_solving.jl b/test/lowering_solving.jl
index 6527edde4b..bfaeee60cf 100644
--- a/test/lowering_solving.jl
+++ b/test/lowering_solving.jl
@@ -8,7 +8,7 @@ eqs = [D(D(x)) ~ σ * (y - x),
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z]
 
-@named sys′ = ODESystem(eqs)
+@named sys′ = ODESystem(eqs, t)
 sys = ode_order_lowering(sys′)
 
 eqs2 = [0 ~ x * y - k,
@@ -46,8 +46,8 @@ eqs = [D(x) ~ σ * (y - x),
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z]
 
-lorenz1 = ODESystem(eqs, name = :lorenz1)
-lorenz2 = ODESystem(eqs, name = :lorenz2)
+lorenz1 = ODESystem(eqs, t, name = :lorenz1)
+lorenz2 = ODESystem(eqs, t, name = :lorenz2)
 
 @variables α(t)
 @parameters γ
diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl
index 3c329c887e..a13625c59a 100644
--- a/test/nonlinearsystem.jl
+++ b/test/nonlinearsystem.jl
@@ -188,7 +188,7 @@ eq = [v1 ~ sin(2pi * t * h)
     v1 - v2 ~ i1
     v2 ~ i2
     i1 ~ i2]
-@named sys = ODESystem(eq)
+@named sys = ODESystem(eq, t)
 @test length(equations(structural_simplify(sys))) == 0
 
 @testset "Issue: 1504" begin
diff --git a/test/odesystem.jl b/test/odesystem.jl
index 65a26f44fd..0a52002f25 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -21,7 +21,7 @@ eqs = [D(x) ~ σ * (y - x),
     D(z) ~ x * y - β * z * κ]
 
 ModelingToolkit.toexpr.(eqs)[1]
-@named de = ODESystem(eqs; defaults = Dict(x => 1))
+@named de = ODESystem(eqs, t; defaults = Dict(x => 1))
 subed = substitute(de, [σ => k])
 ssort(eqs) = sort(eqs, by = string)
 @test isequal(ssort(parameters(subed)), [k, β, ρ])
@@ -29,7 +29,7 @@ ssort(eqs) = sort(eqs, by = string)
     [D(x) ~ k * (y - x)
         D(y) ~ (ρ - z) * x - y
         D(z) ~ x * y - β * κ * z])
-@named des[1:3] = ODESystem(eqs)
+@named des[1:3] = ODESystem(eqs, t)
 @test length(unique(x -> ModelingToolkit.get_tag(x), des)) == 1
 
 @test eval(toexpr(de)) == de
@@ -103,7 +103,7 @@ f = eval(ODEFunctionExpr(de, [x, y, z], [σ, ρ, β], sparsity = false))
 eqs = [D(x) ~ σ * (y - x),
     D(y) ~ x * (ρ - z) - y * t,
     D(z) ~ x * y - β * z * κ]
-@named de = ODESystem(eqs)
+@named de = ODESystem(eqs, t)
 ModelingToolkit.calculate_tgrad(de)
 
 tgrad_oop, tgrad_iip = eval.(ModelingToolkit.generate_tgrad(de))
@@ -119,7 +119,7 @@ tgrad_iip(du, u, p, t)
 eqs = [D(x) ~ σ′ * (y - x),
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z * κ]
-@named de = ODESystem(eqs)
+@named de = ODESystem(eqs, t)
 test_diffeq_inference("global iv-varying", de, t, (x, y, z), (σ′, ρ, β))
 
 f = eval(generate_function(de, [x, y, z], [σ′, ρ, β])[2])
@@ -131,7 +131,7 @@ f(du, [1.0, 2.0, 3.0], [x -> x + 7, 2, 3], 5.0)
 eqs = [D(x) ~ σ(t - 1) * (y - x),
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z * κ]
-@named de = ODESystem(eqs)
+@named de = ODESystem(eqs, t)
 test_diffeq_inference("single internal iv-varying", de, t, (x, y, z), (σ(t - 1), ρ, β))
 f = eval(generate_function(de, [x, y, z], [σ, ρ, β])[2])
 du = [0.0, 0.0, 0.0]
@@ -139,7 +139,7 @@ f(du, [1.0, 2.0, 3.0], [x -> x + 7, 2, 3], 5.0)
 @test du ≈ [11, -3, -7]
 
 eqs = [D(x) ~ x + 10σ(t - 1) + 100σ(t - 2) + 1000σ(t^2)]
-@named de = ODESystem(eqs)
+@named de = ODESystem(eqs, t)
 test_diffeq_inference("many internal iv-varying", de, t, (x,), (σ(t - 2), σ(t^2), σ(t - 1)))
 f = eval(generate_function(de, [x], [σ])[2])
 du = [0.0]
@@ -152,7 +152,7 @@ D2 = D^2
 @variables u(t) uˍtt(t) uˍt(t) xˍt(t)
 eqs = [D3(u) ~ 2(D2(u)) + D(u) + D(x) + 1
     D2(x) ~ D(x) + 2]
-@named de = ODESystem(eqs)
+@named de = ODESystem(eqs, t)
 de1 = ode_order_lowering(de)
 lowered_eqs = [D(uˍtt) ~ 2uˍtt + uˍt + xˍt + 1
     D(xˍt) ~ xˍt + 2
@@ -165,7 +165,7 @@ lowered_eqs = [D(uˍtt) ~ 2uˍtt + uˍt + xˍt + 1
 # issue #219
 @test all(isequal.([ModelingToolkit.var_from_nested_derivative(eq.lhs)[1]
                     for eq in equations(de1)],
-    unknowns(@named lowered = ODESystem(lowered_eqs))))
+    unknowns(@named lowered = ODESystem(lowered_eqs, t))))
 
 test_diffeq_inference("first-order transform", de1, t, [uˍtt, xˍt, uˍt, u, x], [])
 du = zeros(5)
@@ -178,7 +178,7 @@ a = y - x
 eqs = [D(x) ~ σ * a,
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z * κ]
-@named de = ODESystem(eqs)
+@named de = ODESystem(eqs, t)
 generate_function(de, [x, y, z], [σ, ρ, β])
 jac = calculate_jacobian(de)
 @test ModelingToolkit.jacobian_sparsity(de).colptr == sparse(jac).colptr
@@ -190,7 +190,7 @@ f = ODEFunction(complete(de), [x, y, z], [σ, ρ, β])
 _x = y / C
 eqs = [D(x) ~ -A * x,
     D(y) ~ A * x - B * _x]
-@named de = ODESystem(eqs)
+@named de = ODESystem(eqs, t)
 @test begin
     local f, du
     f = eval(generate_function(de, [x, y], [A, B, C])[2])
@@ -231,7 +231,7 @@ ODEFunction(de)(similar(prob.u0), prob.u0, prob.p, 0.1)
 eqs = [D(y₁) ~ -k₁ * y₁ + k₃ * y₂ * y₃,
     0 ~ y₁ + y₂ + y₃ - 1,
     D(y₂) ~ k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃ * κ]
-@named sys = ODESystem(eqs, defaults = [k₁ => 100, k₂ => 3e7, y₁ => 1.0])
+@named sys = ODESystem(eqs, t, defaults = [k₁ => 100, k₂ => 3e7, y₁ => 1.0])
 sys = complete(sys)
 u0 = Pair[]
 push!(u0, y₂ => 0.0)
@@ -272,7 +272,7 @@ sol_dpmap = solve(prob_dpmap, Rodas5())
     function makesys(name)
         @parameters a=1.0
         @variables x(t) = 0.0
-        ODESystem([D(x) ~ -a * x]; name)
+        ODESystem([D(x) ~ -a * x], t; name)
     end
 
     function makecombinedsys()
@@ -321,7 +321,7 @@ end
 eqs = [D(x) ~ σ * (y - x),
     D(y) ~ x - β * y,
     x + z ~ y]
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 @test all(isequal.(unknowns(sys), [x, y, z]))
 @test all(isequal.(parameters(sys), [σ, β]))
 @test equations(sys) == eqs
@@ -331,7 +331,7 @@ eqs = [D(x) ~ σ * (y - x),
 using ModelingToolkit
 @parameters a
 @variables x(t)
-@named sys = ODESystem([D(x) ~ a])
+@named sys = ODESystem([D(x) ~ a], t)
 @test issym(equations(sys)[1].rhs)
 
 # issue 708
@@ -375,7 +375,7 @@ eqs = [
 @parameters r
 @variables x(t)
 eq = D(x) ~ r * x
-@named ode = ODESystem(eq)
+@named ode = ODESystem(eq, t)
 @test equations(ode) == [eq]
 # issue #808
 @testset "Combined system name collisions" begin
@@ -383,7 +383,7 @@ eq = D(x) ~ r * x
         @parameters a
         @variables x(t) f(t)
 
-        ODESystem([D(x) ~ -a * x + f]; name)
+        ODESystem([D(x) ~ -a * x + f], t; name)
     end
 
     function issue808()
@@ -458,7 +458,7 @@ end
 # issue 1109
 let
     @variables x(t)[1:3, 1:3]
-    @named sys = ODESystem(D.(x) .~ x)
+    @named sys = ODESystem(D.(x) .~ x, t)
     @test_nowarn structural_simplify(sys)
 end
 
@@ -493,7 +493,7 @@ function submodel(; name)
     @variables y(t)
     @parameters A[1:5]
     A = collect(A)
-    ODESystem(D(y) ~ sum(A) * y; name = name)
+    ODESystem(D(y) ~ sum(A) * y, t; name = name)
 end
 
 # Build system
@@ -539,7 +539,7 @@ prob = ODEProblem(outersys, [sys.x => 1.0; collect(sys.ms) .=> 1:3], (0, 1.0))
 # observed variable handling
 @variables x(t) RHS(t)
 @parameters τ
-@named fol = ODESystem([D(x) ~ (1 - x) / τ]; observed = [RHS ~ (1 - x) / τ])
+@named fol = ODESystem([D(x) ~ (1 - x) / τ], t; observed = [RHS ~ (1 - x) / τ])
 @test isequal(RHS, @nonamespace fol.RHS)
 RHS2 = RHS
 @unpack RHS = fol
@@ -699,12 +699,6 @@ let
     # @test eltype(prob.p) === Union{Float64, Int}
 end
 
-let
-    @variables t s(t) I(t) r(t)
-    @parameters N
-    @test_throws Any @named tmp = ODESystem([s + I + r ~ N])
-end
-
 let
     @parameters C L R
     @variables q(t) p(t) F(t)
@@ -722,37 +716,35 @@ let
     eqs_to_lhs(eqs) = eq_to_lhs.(eqs)
 
     @parameters σ=10 ρ=28 β=8 / 3 sigma rho beta
-    @variables t2 x(t)=1 y(t)=0 z(t)=0 x2(t2)=1 y2(t2)=0 z2(t2)=0 u(t2)[1:3]
-
-    D2 = Differential(t2)
+    @variables x(t)=1 y(t)=0 z(t)=0 x2(t)=1 y2(t)=0 z2(t)=0 u(t)[1:3]
 
     eqs = [D(x) ~ σ * (y - x),
         D(y) ~ x * (ρ - z) - y,
         D(z) ~ x * y - β * z]
 
     eqs2 = [
-        D2(y2) ~ x2 * (rho - z2) - y2,
-        D2(x2) ~ sigma * (y2 - x2),
-        D2(z2) ~ x2 * y2 - beta * z2,
+        D(y2) ~ x2 * (rho - z2) - y2,
+        D(x2) ~ sigma * (y2 - x2),
+        D(z2) ~ x2 * y2 - beta * z2,
     ]
 
     # array u
-    eqs3 = [D2(u[1]) ~ sigma * (u[2] - u[1]),
-        D2(u[2]) ~ u[1] * (rho - u[3]) - u[2],
-        D2(u[3]) ~ u[1] * u[2] - beta * u[3]]
+    eqs3 = [D(u[1]) ~ sigma * (u[2] - u[1]),
+        D(u[2]) ~ u[1] * (rho - u[3]) - u[2],
+        D(u[3]) ~ u[1] * u[2] - beta * u[3]]
     eqs3 = eqs_to_lhs(eqs3)
 
     eqs4 = [
-        D2(y2) ~ x2 * (rho - z2) - y2,
-        D2(x2) ~ sigma * (y2 - x2),
-        D2(z2) ~ y2 - beta * z2, # missing x2 term
+        D(y2) ~ x2 * (rho - z2) - y2,
+        D(x2) ~ sigma * (y2 - x2),
+        D(z2) ~ y2 - beta * z2, # missing x2 term
     ]
 
-    @named sys1 = ODESystem(eqs)
-    @named sys2 = ODESystem(eqs2)
-    @named sys3 = ODESystem(eqs3, t2)
+    @named sys1 = ODESystem(eqs, t)
+    @named sys2 = ODESystem(eqs2, t)
+    @named sys3 = ODESystem(eqs3, t)
     ssys3 = structural_simplify(sys3)
-    @named sys4 = ODESystem(eqs4)
+    @named sys4 = ODESystem(eqs4, t)
 
     @test ModelingToolkit.isisomorphic(sys1, sys2)
     @test !ModelingToolkit.isisomorphic(sys1, sys3)
@@ -863,7 +855,7 @@ let
         @variables u(t) x(t) v(t)
 
         eqs = [u ~ kx * x + kv * v]
-        ODESystem(eqs; name)
+        ODESystem(eqs, t; name)
     end
 
     @named ctrl = pd_ctrl()
@@ -874,7 +866,7 @@ let
         @variables u(t) x(t) v(t)
 
         eqs = [D(x) ~ v, D(v) ~ u]
-        ODESystem(eqs; name)
+        ODESystem(eqs, t; name)
     end
 
     @named sys = double_int()
@@ -883,7 +875,7 @@ let
 
     connections = [sys.u ~ ctrl.u, ctrl.x ~ sys.x, ctrl.v ~ sys.v]
 
-    @named connected = ODESystem(connections)
+    @named connected = ODESystem(connections, t)
     @named sys_con = compose(connected, sys, ctrl)
 
     sys_simp = structural_simplify(sys_con)
@@ -909,7 +901,7 @@ let
     ∂t = D
     eqs = [∂t(Q) ~ 0.2P
         ∂t(P) ~ -80.0sin(Q)]
-    @test_throws ArgumentError @named sys = ODESystem(eqs)
+    @test_throws ArgumentError @named sys = ODESystem(eqs, t)
 end
 
 @parameters C L R
diff --git a/test/precompile_test/ODEPrecompileTest.jl b/test/precompile_test/ODEPrecompileTest.jl
index caf411f88d..3e25fa21e2 100644
--- a/test/precompile_test/ODEPrecompileTest.jl
+++ b/test/precompile_test/ODEPrecompileTest.jl
@@ -12,7 +12,7 @@ function system(; kwargs...)
         D(y) ~ x * (ρ - z) - y,
         D(z) ~ x * y - β * z]
 
-    @named de = ODESystem(eqs)
+    @named de = ODESystem(eqs, t)
     de = complete(de)
     return ODEFunction(de, [x, y, z], [σ, ρ, β]; kwargs...)
 end
diff --git a/test/reduction.jl b/test/reduction.jl
index 92a73b0102..084b2cdd5e 100644
--- a/test/reduction.jl
+++ b/test/reduction.jl
@@ -240,7 +240,7 @@ vars = @variables x(t) y(t) z(t)
 eqs = [D(x) ~ x
     D(y) ~ y
     D(z) ~ t]
-@named model = ODESystem(eqs)
+@named model = ODESystem(eqs, t)
 sys = structural_simplify(model)
 Js = ModelingToolkit.jacobian_sparsity(sys)
 @test size(Js) == (3, 3)
@@ -257,17 +257,17 @@ ss = alias_elimination(sys)
 
 @variables x(t) y(t)
 @named sys = ODESystem([D(x) ~ 1 - x,
-    D(y) + D(x) ~ 0])
+    D(y) + D(x) ~ 0], t)
 new_sys = alias_elimination(sys)
 @test isempty(observed(new_sys))
 
 @named sys = ODESystem([D(x) ~ x,
-    D(y) + D(x) ~ 0])
+    D(y) + D(x) ~ 0], t)
 new_sys = alias_elimination(sys)
 @test isempty(observed(new_sys))
 
 @named sys = ODESystem([D(x) ~ 1 - x,
-    y + D(x) ~ 0])
+    y + D(x) ~ 0], t)
 new_sys = alias_elimination(sys)
 @test isempty(observed(new_sys))
 
diff --git a/test/serialization.jl b/test/serialization.jl
index 62affc7c67..fb6525b9f7 100644
--- a/test/serialization.jl
+++ b/test/serialization.jl
@@ -3,7 +3,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 
 @variables x(t)
 
-@named sys = ODESystem([D(x) ~ -0.5 * x], defaults = Dict(x => 1.0))
+@named sys = ODESystem([D(x) ~ -0.5 * x], t, defaults = Dict(x => 1.0))
 sys = complete(sys)
 for prob in [
     eval(ModelingToolkit.ODEProblem{false}(sys, nothing, nothing,
diff --git a/test/static_arrays.jl b/test/static_arrays.jl
index 871969c72a..621e15591e 100644
--- a/test/static_arrays.jl
+++ b/test/static_arrays.jl
@@ -8,7 +8,7 @@ eqs = [D(D(x)) ~ σ * (y - x),
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z]
 
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 sys = structural_simplify(sys)
 
 u0 = @SVector [D(x) => 2.0,
diff --git a/test/steadystatesystems.jl b/test/steadystatesystems.jl
index 49d6747f88..d29c809325 100644
--- a/test/steadystatesystems.jl
+++ b/test/steadystatesystems.jl
@@ -6,7 +6,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 @parameters r
 @variables x(t)
 eqs = [D(x) ~ x^2 - r]
-@named de = ODESystem(eqs)
+@named de = ODESystem(eqs, t)
 de = complete(de)
 
 for factor in [1e-1, 1e0, 1e10],
diff --git a/test/structural_transformation/tearing.jl b/test/structural_transformation/tearing.jl
index 8ebb415174..179667e528 100644
--- a/test/structural_transformation/tearing.jl
+++ b/test/structural_transformation/tearing.jl
@@ -152,7 +152,7 @@ newdaesys = structural_simplify(daesys)
 @test equations(newdaesys) == [D(x) ~ z; 0 ~ y + sin(z) - p * t]
 @test equations(tearing_substitution(newdaesys)) == [D(x) ~ z; 0 ~ x + sin(z) - p * t]
 @test isequal(unknowns(newdaesys), [x, z])
-@test isequal(states(newdaesys), [x, z])
+@test isequal(unknowns(newdaesys), [x, z])
 @test_deprecated ODAEProblem(newdaesys, [x => 1.0, z => -0.5π], (0, 1.0), [p => 0.2])
 prob = ODEProblem(newdaesys, [x => 1.0, z => -0.5π], (0, 1.0), [p => 0.2])
 du = [0.0, 0.0];
diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl
index 66c52bc7d8..4cc0d93a4f 100644
--- a/test/symbolic_events.jl
+++ b/test/symbolic_events.jl
@@ -116,14 +116,14 @@ end
 
 ##
 
-@named sys = ODESystem(eqs, continuous_events = [x ~ 1])
+@named sys = ODESystem(eqs, t, continuous_events = [x ~ 1])
 @test getfield(sys, :continuous_events)[] ==
       SymbolicContinuousCallback(Equation[x ~ 1], NULL_AFFECT)
 @test isequal(equations(getfield(sys, :continuous_events))[], x ~ 1)
 fsys = flatten(sys)
 @test isequal(equations(getfield(fsys, :continuous_events))[], x ~ 1)
 
-@named sys2 = ODESystem([D(x) ~ 1], continuous_events = [x ~ 2], systems = [sys])
+@named sys2 = ODESystem([D(x) ~ 1], t, continuous_events = [x ~ 2], systems = [sys])
 @test getfield(sys2, :continuous_events)[] ==
       SymbolicContinuousCallback(Equation[x ~ 2], NULL_AFFECT)
 @test all(ModelingToolkit.continuous_events(sys2) .== [
@@ -191,7 +191,7 @@ sol = solve(prob, Tsit5())
 @test minimum(t -> abs(t - 1), sol.t) < 1e-10 # test that the solver stepped at the first root
 @test minimum(t -> abs(t - 2), sol.t) < 1e-10 # test that the solver stepped at the second root
 
-@named sys = ODESystem(eqs, continuous_events = [x ~ 1, x ~ 2]) # two root eqs using the same unknown
+@named sys = ODESystem(eqs, t, continuous_events = [x ~ 1, x ~ 2]) # two root eqs using the same unknown
 sys = complete(sys)
 prob = ODEProblem(sys, Pair[], (0.0, 3.0))
 @test get_callback(prob) isa ModelingToolkit.DiffEqCallbacks.VectorContinuousCallback
@@ -288,7 +288,7 @@ eq = [vs ~ sin(2pi * t)
     D(v) ~ vs - v
     D(vmeasured) ~ 0.0]
 ev = [sin(20pi * t) ~ 0.0] => [vmeasured ~ v]
-@named sys = ODESystem(eq, continuous_events = ev)
+@named sys = ODESystem(eq, t, continuous_events = ev)
 sys = structural_simplify(sys)
 prob = ODEProblem(sys, zeros(2), (0.0, 5.1))
 sol = solve(prob, Tsit5())
diff --git a/test/units.jl b/test/units.jl
index d1a6c57640..bdaa495f5a 100644
--- a/test/units.jl
+++ b/test/units.jl
@@ -37,13 +37,13 @@ D = Differential(t)
 @test UMT.get_unit(1.0^(t / τ)) == UMT.unitless
 @test UMT.get_unit(exp(t / τ)) == UMT.unitless
 @test UMT.get_unit(sin(t / τ)) == UMT.unitless
-@test UMT.get_unit(sin(1u"rad")) == UMT.unitless
+@test UMT.get_unit(sin(1 * u"rad")) == UMT.unitless
 @test UMT.get_unit(t^2) == u"ms^2"
 
 eqs = [D(E) ~ P - E / τ
     0 ~ P]
 @test UMT.validate(eqs)
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 
 @test !UMT.validate(D(D(E)) ~ P)
 @test !UMT.validate(0 ~ P + E * τ)
@@ -53,13 +53,13 @@ eqs = [D(E) ~ P - E / τ
 ODESystem(eqs, t, [E, P, t], [τ], name = :sys, checks = MT.CheckUnits)
 eqs = [D(E) ~ P - E / τ
     0 ~ P + E * τ]
-@test_throws MT.ValidationError ODESystem(eqs, name = :sys, checks = MT.CheckAll)
-@test_throws MT.ValidationError ODESystem(eqs, name = :sys, checks = true)
-ODESystem(eqs, name = :sys, checks = MT.CheckNone)
-ODESystem(eqs, name = :sys, checks = false)
-@test_throws MT.ValidationError ODESystem(eqs, name = :sys,
+@test_throws MT.ValidationError ODESystem(eqs, t, name = :sys, checks = MT.CheckAll)
+@test_throws MT.ValidationError ODESystem(eqs, t, name = :sys, checks = true)
+ODESystem(eqs, t, name = :sys, checks = MT.CheckNone)
+ODESystem(eqs, t, name = :sys, checks = false)
+@test_throws MT.ValidationError ODESystem(eqs, t, name = :sys,
     checks = MT.CheckComponents | MT.CheckUnits)
-@named sys = ODESystem(eqs, checks = MT.CheckComponents)
+@named sys = ODESystem(eqs, t, checks = MT.CheckComponents)
 @test_throws MT.ValidationError ODESystem(eqs, t, [E, P, t], [τ], name = :sys,
     checks = MT.CheckUnits)
 
@@ -98,7 +98,7 @@ bad_length_eqs = [connect(op, lp)]
 @parameters v[1:3]=[1, 2, 3] [unit = u"m/s"]
 D = Differential(t)
 eqs = D.(x) .~ v
-ODESystem(eqs, name = :sys)
+ODESystem(eqs, t, name = :sys)
 
 # Nonlinear system
 @parameters a [unit = u"kg"^-1]
@@ -135,12 +135,12 @@ noiseeqs = [0.1u"MW" 0.1u"MW"
 D = Differential(t)
 eqs = [D(L) ~ v,
     V ~ L^3]
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 sys_simple = structural_simplify(sys)
 
 eqs = [D(V) ~ r,
     V ~ L^3]
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 sys_simple = structural_simplify(sys)
 
 @variables V [unit = u"m"^3] L [unit = u"m"]