diff --git a/Project.toml b/Project.toml
index 752d22ae6..630b3d6d5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -52,7 +52,7 @@ Random = "1.6"
 RandomNumbers = "1.5.3"
 RecursiveArrayTools = "2, 3"
 Reexport = "0.2, 1.0"
-SciMLBase = "2.51"
+SciMLBase = "2.59.2"
 SciMLOperators = "0.2.9, 0.3"
 SparseArrays = "1.6"
 SparseDiffTools = "2"
diff --git a/src/solve.jl b/src/solve.jl
index 28911c460..d7bc2f8c8 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -230,6 +230,7 @@ function DiffEqBase.__init(
   d_discontinuities_internal = OrdinaryDiffEq.initialize_d_discontinuities(tType, d_discontinuities, tspan)
 
   ### Algorithm-specific defaults ###
+  save_idxs, saved_subsystem = SciMLBase.get_save_idxs_and_saved_subsystem(prob, save_idxs)
   # if save_idxs === nothing
   #   ksEltype = Vector{rateType}
   # else
@@ -538,12 +539,12 @@ function DiffEqBase.__init(
   if alg isa Union{StochasticDiffEqCompositeAlgorithm,
                           StochasticDiffEqRODECompositeAlgorithm}
     sol = DiffEqBase.build_solution(prob,alg,ts,timeseries,W=W,
-                                    stats = stats,
+                                    stats = stats, saved_subsystem = saved_subsystem,
                                     calculate_error = false, alg_choice=alg_choice,
                                     interp = id, dense = dense, seed = _seed)
   else
     sol = DiffEqBase.build_solution(prob,alg,ts,timeseries,W=W,
-                                    stats = stats,
+                                    stats = stats, saved_subsystem = saved_subsystem,
                                     calculate_error = false,
                                     interp = id, dense = dense, seed = _seed)
   end