feat: allow interpolating discretes past the final time point

Author: AayushSabharwal

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

test/downstream/solution_interface.jl

@@ -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