Performant way to remake ODEProblem with updated external input function? about modelingtoolkit.jl HOT 2 CLOSED

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using DifferentialEquations
using DataInterpolations
using ForwardDiff

P0 = 2.0 # test value for the parameter P

# F(t) is an external function that depends on the parameter P
function create_F(P)
    ts = range(0.0, 1.0, step=0.1)
    Fs = ts .^ P # toy example: F = t^P
    return CubicSpline(Fs, ts)
end
F(t, spline) = spline(t) # intermediate function for F(t) (possible to avoid this?)

# an ODE that depends on the parameter P both explicitly, and implicitly through F(t)
@parameters P F_spline
@variables x(t)
@register_symbolic F(t, spline) # following https://docs.sciml.ai/ModelingToolkit/stable/tutorials/ode_modeling/#Specifying-a-time-variable-forcing-function
sys = structural_simplify(ODESystem([D(x) ~ P * F(t, F_spline)], t, [x], [P, F_spline]; name = :sys))
prob = ODEProblem(sys, [x => 0.0], (0.0, 1.0), [P => NaN, F_spline => NaN]) # uninitialized problem

# solution of the ODE for a given value of the parameter P (must be fast!)
function solve_instance_fast(P)
    F_instance = create_F(P)
    prob_instance = remake(prob; p = [sys.P => P, sys.F_spline => F_instance]) # update P and spline for F(t)
    return solve(prob_instance)
end

x_at_1(P) = solve_instance_fast(P)(1.0; idxs=x) # analytical result: x(1) = P/(P+1)
ForwardDiff.derivative(x_at_1, P0) # analytical result: d(x(t=1))/dP = 1/(P+1)^2