diff --git a/docs/Project.toml b/docs/Project.toml
index fa8663361..9d2ea8cf8 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -62,7 +62,7 @@ ModelingToolkit = "8"
 ODEProblemLibrary = "0.1"
 Optimization = "3"
 OptimizationNLopt = "0.1, 0.2"
-OrdinaryDiffEq = "6.31"
+OrdinaryDiffEq = "6.76"
 Plots = "1"
 RecursiveArrayTools = "2, 3"
 SDEProblemLibrary = "0.1"
diff --git a/docs/src/examples/kepler_problem.md b/docs/src/examples/kepler_problem.md
index 053c0a559..0a6008b40 100644
--- a/docs/src/examples/kepler_problem.md
+++ b/docs/src/examples/kepler_problem.md
@@ -10,8 +10,7 @@ Also, we know that
 $${\displaystyle {\frac {\mathrm {d} {\boldsymbol {p}}}{\mathrm {d} t}}=-{\frac {\partial {\mathcal {H}}}{\partial {\boldsymbol {q}}}}\quad ,\quad {\frac {\mathrm {d} {\boldsymbol {q}}}{\mathrm {d} t}}=+{\frac {\partial {\mathcal {H}}}{\partial {\boldsymbol {p}}}}}$$
 
 ```@example kepler
-using OrdinaryDiffEq, LinearAlgebra, ForwardDiff, Plots;
-gr();
+using OrdinaryDiffEq, LinearAlgebra, ForwardDiff, Plots
 H(q, p) = norm(p)^2 / 2 - inv(norm(q))
 L(q, p) = q[1] * p[2] - p[1] * q[2]
 
@@ -30,12 +29,12 @@ sol = solve(prob, KahanLi6(), dt = 1 // 10);
 Let's plot the orbit and check the energy and angular momentum variation. We know that energy and angular momentum should be constant, and they are also called first integrals.
 
 ```@example kepler
-plot_orbit(sol) = plot(sol, vars = (3, 4), lab = "Orbit", title = "Kepler Problem Solution")
+plot_orbit(sol) = plot(sol, idxs = (3, 4), lab = "Orbit", title = "Kepler Problem Solution")
 
 function plot_first_integrals(sol, H, L)
-    plot(initial_first_integrals[1] .- map(u -> H(u[2, :], u[1, :]), sol.u),
+    plot(initial_first_integrals[1] .- map(u -> H(u.x[2], u.x[1]), sol.u),
         lab = "Energy variation", title = "First Integrals")
-    plot!(initial_first_integrals[2] .- map(u -> L(u[2, :], u[1, :]), sol.u),
+    plot!(initial_first_integrals[2] .- map(u -> L(u.x[2], u.x[1]), sol.u),
         lab = "Angular momentum variation")
 end
 analysis_plot(sol, H, L) = plot(plot_orbit(sol), plot_first_integrals(sol, H, L))