diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index e643be904d..c9ed8e21f1 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -320,8 +320,8 @@ for traitT in [
             elseif is_timeseries_parameter(sys, s)
                 push!(ts_idxs, timeseries_parameter_index(sys, s).timeseries_idx)
             elseif is_time_dependent(sys) && iscall(s) && issym(operation(s)) &&
-                   is_variable(sys, operation(s)(get_iv(sys)))
-                # DDEs case, to detect x(t - k)
+                   is_variable(sys, operation(s)(get_iv(sys), arguments(s)[2:end]...))
+                # DDEs case, to detect x(t - k), and x(t - k, arg2, arg3, ...)
                 push!(ts_idxs, ContinuousTimeseries())
             else
                 if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing
diff --git a/test/dde.jl b/test/dde.jl
index df4acf377a..cf65f10696 100644
--- a/test/dde.jl
+++ b/test/dde.jl
@@ -76,8 +76,9 @@ prob = SDDEProblem(hayes_modelf, hayes_modelg, [1.0], h, tspan, pmul;
     constant_lags = (pmul[1],));
 sol = solve(prob, RKMil(), seed = 100)
 
-@variables x(..) delx(t)
 @parameters a=-4.0 b=-2.0 c=10.0 α=-1.3 β=-1.2 γ=1.1
+@variables x(..)
+@variables delx(t, a, b, c, α, β, γ) # equivalent to just delx(t)
 @brownians η
 τ = 1.0
 eqs = [D(x(t)) ~ a * x(t) + b * x(t - τ) + c + (α * x(t) + γ) * η, delx ~ x(t - τ)]
@@ -88,6 +89,8 @@ eqs = [D(x(t)) ~ a * x(t) + b * x(t - τ) + c + (α * x(t) + γ) * η, delx ~ x(
 @test equations(sys) == [D(x(t)) ~ a * x(t) + b * x(t - τ) + c]
 @test isequal(ModelingToolkit.get_noise_eqs(sys), [α * x(t) + γ;;])
 prob_mtk = SDDEProblem(sys, [x(t) => 1.0 + t], tspan; constant_lags = (τ,));
+@test prob_mtk[delx] isa Float64
+@test prob_mtk[x(t - τ)] isa Float64
 @test_nowarn sol_mtk = solve(prob_mtk, RKMil(), seed = 100)
 
 prob_sa = SDDEProblem(
diff --git a/test/odesystem.jl b/test/odesystem.jl
index 28bc1d7db0..571beb9001 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -1573,3 +1573,13 @@ end
     prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0))
     @test prob.problem_type == "A"
 end
+
+# https://github.com/SciML/ModelingToolkit.jl/issues/3737
+@testset "ODE with multivariate variables" begin
+    @parameters k
+    @variables x(t, k) y(t, k) # equivalent to x(t) and y(t)
+    @mtkbuild sys = System([D(x) ~ 0, y ~ x + 1], t)
+    prob = ODEProblem(sys, [x => 0.0], (0.0, 1.0))
+    @test prob[x] == 0.0 # unknown
+    @test prob[y] == 1.0 # observed
+end