diff --git a/Project.toml b/Project.toml
index 6bca8e7ec..b51c6edec 100644
--- a/Project.toml
+++ b/Project.toml
@@ -97,6 +97,7 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PartialFunctions = "570af359-4316-4cb7-8c74-252c00c2016b"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 RCall = "6f49c342-dc21-5d91-9882-a32aef131414"
@@ -104,7 +105,8 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Pkg", "PyCall", "PythonCall", "SafeTestsets", "Test", "StaticArrays", "StochasticDiffEq", "Aqua", "Zygote", "PartialFunctions", "DataFrames", "ModelingToolkit", "OrdinaryDiffEq"]
+test = ["Pkg", "Plots", "UnicodePlots", "PyCall", "PythonCall", "SafeTestsets", "Test", "StaticArrays", "StochasticDiffEq", "Aqua", "Zygote", "PartialFunctions", "DataFrames", "ModelingToolkit", "OrdinaryDiffEq"]
diff --git a/src/solutions/solution_interface.jl b/src/solutions/solution_interface.jl
index a640dcda7..442e1fe9e 100644
--- a/src/solutions/solution_interface.jl
+++ b/src/solutions/solution_interface.jl
@@ -490,9 +490,10 @@ function solplot_vecs_and_labels(dims, vars, plott, sol, plot_analytic,
                         push!(strs, "u[$(x[j])]")
                     end
                 else
+                    global Main.debug[] = plot_analytic_timeseries
                     _tmp = Vector{eltype(sol[1])}(undef, length(plot_analytic_timeseries))
                     for n in 1:length(plot_analytic_timeseries)
-                        _tmp[n] = plot_analytic_timeseries[n][j]
+                        _tmp[n] = plot_analytic_timeseries[n][x[j]]
                     end
                     push!(tmp, _tmp)
                     if !isempty(varsyms) && x[j] isa Integer
diff --git a/test/downstream/solution_interface.jl b/test/downstream/solution_interface.jl
index da454424a..52ea720f4 100644
--- a/test/downstream/solution_interface.jl
+++ b/test/downstream/solution_interface.jl
@@ -1,5 +1,6 @@
 using ModelingToolkit, OrdinaryDiffEq, RecursiveArrayTools, StochasticDiffEq, Test
 using ModelingToolkit: t_nounits as t, D_nounits as D
+using Plots: Plots, plot
 
 ### Tests on non-layered model (everything should work). ###
 
@@ -25,6 +26,29 @@ sol = solve(oprob, Rodas4())
 @test_throws Exception sol[population_model.a]
 @test_throws Exception sol[:a]
 
+@testset "plot ODE solution" begin
+    Plots.unicodeplots()
+    f = ODEFunction((u, p, t) -> -u, analytic = (u0, p, t) -> u0 * exp(-t))
+
+    # scalar
+    ode = ODEProblem(f, 1.0, (0.0, 1.0))
+    sol = solve(ode, Tsit5())
+    @test_nowarn plot(sol)
+    @test_nowarn plot(sol; plot_analytic = true)
+
+    # vector
+    ode = ODEProblem(f, [1.0, 2.0], (0.0, 1.0))
+    sol = solve(ode, Tsit5())
+    @test_nowarn plot(sol)
+    @test_nowarn plot(sol; plot_analytic = true)
+
+    # matrix
+    ode = ODEProblem(f, [1.0, 2.0; 3.0, 4.0], (0.0, 1.0))
+    sol = solve(ode, Tsit5())
+    @test_nowarn plot(sol)
+    @test_nowarn plot(sol; plot_analytic = true)
+end
+
 # Tests on SDEProblem
 noiseeqs = [0.1 * s1,
     0.1 * s2]