diff --git a/docs/src/examples/higher_order.md b/docs/src/examples/higher_order.md
index ef070db299..d4b2203292 100644
--- a/docs/src/examples/higher_order.md
+++ b/docs/src/examples/higher_order.md
@@ -55,5 +55,5 @@ p  = [σ => 28.0,
 tspan = (0.0,100.0)
 prob = ODEProblem(sys,u0,tspan,p,jac=true)
 sol = solve(prob,Tsit5())
-using Plots; plot(sol,vars=(x,y))
+using Plots; plot(sol,idxs=(x,y))
 ```
diff --git a/docs/src/examples/modelingtoolkitize_index_reduction.md b/docs/src/examples/modelingtoolkitize_index_reduction.md
index af9ce95f44..bcb1611518 100644
--- a/docs/src/examples/modelingtoolkitize_index_reduction.md
+++ b/docs/src/examples/modelingtoolkitize_index_reduction.md
@@ -32,7 +32,7 @@ traced_sys = modelingtoolkitize(pendulum_prob)
 pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
 prob = ODAEProblem(pendulum_sys, [], tspan)
 sol = solve(prob, Tsit5(),abstol=1e-8,reltol=1e-8)
-plot(sol, vars=states(traced_sys))
+plot(sol, idxs=states(traced_sys))
 ```
 
 ## Explanation
@@ -150,7 +150,7 @@ prob = ODEProblem(pendulum_sys, Pair[], tspan)
 sol = solve(prob, Rodas4())
 
 using Plots
-plot(sol, vars=states(traced_sys))
+plot(sol, idxs=states(traced_sys))
 ```
 
 Note that plotting using `states(traced_sys)` is done so that any
@@ -169,7 +169,7 @@ traced_sys = modelingtoolkitize(pendulum_prob)
 pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
 prob = ODAEProblem(pendulum_sys, Pair[], tspan)
 sol = solve(prob, Tsit5(),abstol=1e-8,reltol=1e-8)
-plot(sol, vars=states(traced_sys))
+plot(sol, idxs=states(traced_sys))
 ```
 
 And there you go: this has transformed the model from being too hard to
diff --git a/docs/src/examples/spring_mass.md b/docs/src/examples/spring_mass.md
index 613145c5b8..bdb652bfd9 100644
--- a/docs/src/examples/spring_mass.md
+++ b/docs/src/examples/spring_mass.md
@@ -197,5 +197,5 @@ sol[mass.pos[1]]
 We can also plot the path of the mass:
 
 ```@example component
-plot(sol, vars = (mass.pos[1], mass.pos[2]))
+plot(sol, idxs = (mass.pos[1], mass.pos[2]))
 ```
diff --git a/docs/src/tutorials/acausal_components.md b/docs/src/tutorials/acausal_components.md
index 6331fa503c..6703c55757 100644
--- a/docs/src/tutorials/acausal_components.md
+++ b/docs/src/tutorials/acausal_components.md
@@ -376,5 +376,5 @@ sol[resistor.v]
 or we can plot the timeseries of the resistor's voltage:
 
 ```@example acausal
-plot(sol, vars=[resistor.v])
+plot(sol, idxs=[resistor.v])
 ```
diff --git a/test/lowering_solving.jl b/test/lowering_solving.jl
index 543f5bdab3..2d544dfaf0 100644
--- a/test/lowering_solving.jl
+++ b/test/lowering_solving.jl
@@ -35,7 +35,7 @@ probexpr = ODEProblemExpr(sys, u0, tspan, p, jac = true)
 sol = solve(prob, Tsit5())
 solexpr = solve(eval(prob), Tsit5())
 @test all(x -> x == 0, Array(sol - solexpr))
-#using Plots; plot(sol,vars=(:x,:y))
+#using Plots; plot(sol,idxs=(:x,:y))
 
 @parameters t σ ρ β
 @variables x(t) y(t) z(t)
@@ -73,4 +73,4 @@ tspan = (0.0, 100.0)
 prob = ODEProblem(connected, u0, tspan, p)
 sol = solve(prob, Rodas5())
 @test maximum(sol[2, :] + sol[6, :] + 2sol[1, :]) < 1e-12
-#using Plots; plot(sol,vars=(:α,Symbol(lorenz1.x),Symbol(lorenz2.y)))
+#using Plots; plot(sol,idxs=(:α,Symbol(lorenz1.x),Symbol(lorenz2.y)))
diff --git a/test/structural_transformation/index_reduction.jl b/test/structural_transformation/index_reduction.jl
index f4ebb27d47..44b4b8c534 100644
--- a/test/structural_transformation/index_reduction.jl
+++ b/test/structural_transformation/index_reduction.jl
@@ -70,7 +70,7 @@ prob = ODEProblem(ODEFunction(first_order_idx1_pendulum),
                   [1, 9.8],
                   mass_matrix = calculate_massmatrix(first_order_idx1_pendulum))
 sol = solve(prob, Rodas5());
-#plot(sol, vars=(1, 2))
+#plot(sol, idxs=(1, 2))
 
 new_sys = dae_index_lowering(ModelingToolkit.ode_order_lowering(pendulum2))
 
@@ -83,7 +83,7 @@ prob_auto = ODEProblem(new_sys,
                        (0, 100.0),
                        [1, 9.8])
 sol = solve(prob_auto, Rodas5());
-#plot(sol, vars=(x, y))
+#plot(sol, idxs=(x, y))
 
 # Define some variables
 @parameters t L g
@@ -116,7 +116,7 @@ p = [
 
 prob_auto = ODEProblem(new_sys, u0, (0.0, 10.0), p)
 sol = solve(prob_auto, Rodas5());
-#plot(sol, vars=(D(x), y))
+#plot(sol, idxs=(D(x), y))
 
 let pss_pendulum2 = partial_state_selection(pendulum2)
     @test length(equations(pss_pendulum2)) <= 6