diff --git a/src/remake.jl b/src/remake.jl
index f10e368e55..e7d9903186 100644
--- a/src/remake.jl
+++ b/src/remake.jl
@@ -1083,6 +1083,10 @@ calling `SymbolicIndexingInterface.symbolic_container`, provided for dispatch. R
 the updated `newu0` and `newp`.
 """
 function late_binding_update_u0_p(prob, root_indp, u0, p, t0, newu0, newp)
+    if hasmethod(symbolic_container, Tuple{typeof(root_indp)}) &&
+       (sc = symbolic_container(root_indp)) !== root_indp
+        return late_binding_update_u0_p(prob, sc, u0, p, t0, newu0, newp)
+    end
     return newu0, newp
 end
 
@@ -1094,10 +1098,6 @@ Calls `late_binding_update_u0_p(prob, root_indp, u0, p, t0, newu0, newp)` after
 """
 function late_binding_update_u0_p(prob, u0, p, t0, newu0, newp)
     root_indp = prob
-    while hasmethod(symbolic_container, Tuple{typeof(root_indp)}) &&
-        (sc = symbolic_container(root_indp)) !== root_indp
-        root_indp = sc
-    end
     return late_binding_update_u0_p(prob, root_indp, u0, p, t0, newu0, newp)
 end
 
diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl
index 83b9580e49..46ed79bae2 100644
--- a/src/scimlfunctions.jl
+++ b/src/scimlfunctions.jl
@@ -2808,7 +2808,8 @@ function unwrapped_f(f::NonlinearFunction, newf = unwrapped_f(f.f))
     end
 end
 
-@add_kwonly function SplitFunction(f1, f2, mass_matrix, _func_cache, analytic, tgrad, jac, jvp,
+@add_kwonly function SplitFunction(
+        f1, f2, mass_matrix, _func_cache, analytic, tgrad, jac, jvp,
         vjp, jac_prototype, W_prototype, sparsity, Wfact, Wfact_t, paramjac,
         observed, colorvec, sys, initializeprob = nothing, update_initializeprob! = nothing,
         initializeprobmap = nothing, initializeprobpmap = nothing, initialization_data = nothing, nlprob_data = nothing)
@@ -3259,7 +3260,8 @@ function SDEFunction(f, g; kwargs...)
 end
 SDEFunction(f::SDEFunction; kwargs...) = f
 
-@add_kwonly function SplitSDEFunction(f1, f2, g, mass_matrix, _func_cache, analytic, tgrad, jac,
+@add_kwonly function SplitSDEFunction(
+        f1, f2, g, mass_matrix, _func_cache, analytic, tgrad, jac,
         jvp, vjp,
         jac_prototype, Wfact, Wfact_t, paramjac, observed,
         colorvec, sys, initialization_data = nothing)
@@ -3344,7 +3346,8 @@ function SplitSDEFunction{iip}(f1, f2, g; kwargs...) where {iip}
 end
 SplitSDEFunction(f::SplitSDEFunction; kwargs...) = f
 
-@add_kwonly function DynamicalSDEFunction(f1, f2, g, mass_matrix, _func_cache, analytic, tgrad,
+@add_kwonly function DynamicalSDEFunction(
+        f1, f2, g, mass_matrix, _func_cache, analytic, tgrad,
         jac, jvp, vjp,
         jac_prototype, Wfact, Wfact_t, paramjac,
         observed, colorvec,
diff --git a/src/solutions/ode_solutions.jl b/src/solutions/ode_solutions.jl
index 149cc18e33..dcae2c08f8 100644
--- a/src/solutions/ode_solutions.jl
+++ b/src/solutions/ode_solutions.jl
@@ -287,8 +287,7 @@ function (sol::AbstractODESolution)(t::Number, ::Type{deriv}, idxs,
         ps = parameter_values(discs)
         for ts_idx in eachindex(discs)
             partition = discs[ts_idx]
-            interp_val = ConstantInterpolation(partition.t, partition.u)(
-                t, nothing, deriv, nothing, continuity)
+            interp_val = _hold_discrete(partition.u, partition.t, t)
             ps = with_updated_parameter_timeseries_values(sol, ps, ts_idx => interp_val)
         end
     end
@@ -312,8 +311,7 @@ function (sol::AbstractODESolution)(t::Number, ::Type{deriv}, idxs::AbstractVect
         ps = parameter_values(discs)
         for ts_idx in eachindex(discs)
             partition = discs[ts_idx]
-            interp_val = ConstantInterpolation(partition.t, partition.u)(
-                t, nothing, deriv, nothing, continuity)
+            interp_val = _hold_discrete(partition.u, partition.t, t)
             ps = with_updated_parameter_timeseries_values(sol, ps, ts_idx => interp_val)
         end
     end
diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml
index d148b10b7b..c390a258f4 100644
--- a/test/downstream/Project.toml
+++ b/test/downstream/Project.toml
@@ -34,7 +34,7 @@ DelayDiffEq = "5"
 DiffEqCallbacks = "3, 4"
 ForwardDiff = "0.10"
 JumpProcesses = "9.10"
-ModelingToolkit = "9.56"
+ModelingToolkit = "9.64.1"
 ModelingToolkitStandardLibrary = "2.7"
 NonlinearSolve = "2, 3, 4"
 Optimization = "4"
diff --git a/test/downstream/solution_interface.jl b/test/downstream/solution_interface.jl
index 0235cc187d..6d363967d1 100644
--- a/test/downstream/solution_interface.jl
+++ b/test/downstream/solution_interface.jl
@@ -215,7 +215,7 @@ end
             @test sol[x] == xvals
             @test is_parameter(sol, p)
             @test parameter_index(sol, p) == parameter_index(sys, p)
-            @test isequal(only(parameter_symbols(sol)), p)
+            @test any(isequal(p), parameter_symbols(sol))
             @test is_independent_variable(sol, t)
 
             tmp = copy(prob.u0)
@@ -341,3 +341,12 @@ end
         @test _ss isa SciMLBase.SavedSubsystem
     end
 end
+
+@testset "Interpolation after final discrete save" begin
+    @variables x(t) y(t)
+    @parameters start
+    @mtkbuild sys=ODESystem([D(x) ~ y, y ~ ifelse(t < start, 1.0, 2.0)], t) additional_passes=[ModelingToolkit.IfLifting]
+    prob = ODEProblem(sys, [x => 0.0], (0.0, 1.0), [start => 0.5])
+    sol = solve(prob)
+    @test sol(0.6, idxs = y) ≈ 2.0
+end