Problems with solving one ODE with the solution of another

[hersle]

Here is a simple ODE system (from some perturbation theory toy problem; the real problem is bigger):

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

@variables y(t)[0:1]
eqs0 = [D(y[0]) ~ -1] # 0th order equations
eqs1 = [D(y[1]) ~ 2*y[0]] # 1st order equations
ics0 = [y[0] => 0] # 0th order initial conditions
ics1 = [y[1] => 0] # 1st order initial conditions

# 0th+1st order system
@named sys = ODESystem([eqs0; eqs1], t; defaults = [ics0; ics1])

I can easily solve the whole system "simultaneously":

ssys = structural_simplify(sys)
prob = ODEProblem(ssys, [], (0, 3))
sol = solve(prob)

---

But now I want to exploit this property of the system: y[0] is independent of y[1] (but y[1] depends on y[0]). This always happens in perturbation theory. I want to solve for only y[0] first, and then use that ODESolution to solve for only y[1]. I split the "simultaneous" system into two "sequential" ones (my dream is to make a function that does this for a general system, i.e. a block lower triangular transformation):

@register_symbolic evalsol(sol::ODESolution, t::AbstractFloat)
evalsol(sol::ODESolution, t::AbstractFloat) = sol(t, idxs=y[0])

# 0th order system
@named sys0 = ODESystem(eqs0, t; defaults = ics0)

# 1st order system: replace 0th order ODE by its solution
@parameters sol0::ODESolution
@named sys1 = ODESystem([y[0] ~ evalsol(sol0, t); eqs1], t; defaults = ics1)

I now solve them "sequentially", and everything works!

ssys0, ssys1 = structural_simplify(sys0), structural_simplify(sys1)

prob0 = ODEProblem(ssys0, [], (0, 3))
sol0 = solve(prob0)

prob1 = ODEProblem(ssys1, [], (0, 3), [ssys1.sol0 => sol0])
sol1 = solve(prob1)

using Test
@test sol(sol.t, idxs=y[0]) ≈ sol0(sol.t, idxs=y[0]) # passes!
@test sol(sol.t, idxs=y[0]) ≈ sol1(sol.t, idxs=y[0]) # passes!
@test sol(sol.t, idxs=y[1]) ≈ sol1(sol.t, idxs=y[1]) # passes!

But there are some problems that block me from doing this generally (next post).

[hersle]

The first problem is that I want to get rid of the evalsol() workaround that hardcodes idxs=y[0]. Ideally, I'd like to use callable parameters, but to start with I try

@register_symbolic evalsol(sol::ODESolution, t::AbstractFloat, idx::Num)
evalsol(sol::ODESolution, t::AbstractFloat, idx::Num) = sol(t, idxs=idx)

# 1st order system: replace 0th order ODE by its solution
@parameters sol0::ODESolution
@named sys1 = ODESystem([y[0] ~ evalsol(sol0, t, y[0]); eqs1], t; defaults = ics1)
ssys1 = structural_simplify(sys1)

Now the problem is that the 0th order equation in ssys1 is no longer an observed equation, which it should be:

equations(ssys1)

2-element Vector{Equation}:
 Differential(t)((y(t))[1]) ~ 2(y(t))[0]
 0 ~ -(y(t))[0] + evalsol(sol0, t, (y(t))[0]) # <-- y[0] should be observed!

That makes ssys1 an unbalanced/unsolvable system (in the same way as above). Is there a way around this?