diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 958000ec3e..d91659c0d7 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -2878,6 +2878,15 @@ function Base.eltype(::Type{<:TreeIterator{ModelingToolkit.AbstractSystem}})
     ModelingToolkit.AbstractSystem
 end
 
+function check_array_equations_unknowns(eqs, dvs)
+    if any(eq -> Symbolics.isarraysymbolic(eq.lhs), eqs)
+        throw(ArgumentError("The system has array equations. Call `structural_simplify` to handle such equations or scalarize them manually."))
+    end
+    if any(x -> Symbolics.isarraysymbolic(x), dvs)
+        throw(ArgumentError("The system has array unknowns. Call `structural_simplify` to handle this or scalarize them manually."))
+    end
+end
+
 function check_eqs_u0(eqs, dvs, u0; check_length = true, kwargs...)
     if u0 !== nothing
         if check_length
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 06feaac688..a7ae411def 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -815,6 +815,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
     ps = parameters(sys)
     iv = get_iv(sys)
 
+    check_array_equations_unknowns(eqs, dvs)
     # TODO: Pass already computed information to varmap_to_vars call
     # in process_u0? That would just be a small optimization
     varmap = u0map === nothing || isempty(u0map) || eltype(u0map) <: Number ?
diff --git a/test/odesystem.jl b/test/odesystem.jl
index 47c8ad9661..bdaf489a32 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -1404,3 +1404,34 @@ end
     @test !ModelingToolkit.is_dde(sys)
     @test is_markovian(sys)
 end
+
+@testset "Issue #2597" begin
+    @variables x(t)[1:2]=ones(2) y(t)=1.0
+
+    for eqs in [D(x) ~ x, collect(D(x) .~ x)]
+        for dvs in [[x], collect(x)]
+            @named sys = ODESystem(eqs, t, dvs, [])
+            sys = complete(sys)
+            if eqs isa Vector && length(eqs) == 2 && length(dvs) == 2
+                @test_nowarn ODEProblem(sys, [], (0.0, 1.0))
+            else
+                @test_throws [
+                    r"array (equations|unknowns)", "structural_simplify", "scalarize"] ODEProblem(
+                    sys, [], (0.0, 1.0))
+            end
+        end
+    end
+    for eqs in [[D(x) ~ x, D(y) ~ y], [collect(D(x) .~ x); D(y) ~ y]]
+        for dvs in [[x, y], [x..., y]]
+            @named sys = ODESystem(eqs, t, dvs, [])
+            sys = complete(sys)
+            if eqs isa Vector && length(eqs) == 3 && length(dvs) == 3
+                @test_nowarn ODEProblem(sys, [], (0.0, 1.0))
+            else
+                @test_throws [
+                    r"array (equations|unknowns)", "structural_simplify", "scalarize"] ODEProblem(
+                    sys, [], (0.0, 1.0))
+            end
+        end
+    end
+end