Niche case of `plotsol; idxs = [...])` fails

[TorkelE]

I think this might be a case of regression. I this used to work in Catalyst tests, but stopped working at some point and I didn't have time to figure it out then. Going through broken Catalyst tests and noted this one had new issue.

MWE:

using ModelingToolkit, OrdinaryDiffEq, Plots
@variables t X1(t) X2(t) Y1(t) Y2(t)
@parameters p d
D = Differential(t)
eqs = [
    D(X1) ~ p - d*X1,
    D(Y1) ~ p - d*Y1,
    X2 ~ X1 + 1,
    Y2 ~ Y1 + 1
]
@mtkbuild osys = ODESystem(eqs, t)

u0 = [X1 => 2.0, Y1 => 1.0]
ps = [p => 1.0, d => 0.2]
oprob = ODEProblem(osys, u0, 1.0, ps)
sol = solve(oprob)
plot(sol; idxs = [:X2]) # Works
plot(sol; idxs = [:Y2]) # Works
plot(sol; idxs = [X2, Y2]) # Works.
plot(sol; idxs = [osys.X2, osys.Y2]) # Works.
plot(sol; idxs = [:X2, :Y2]) # `ERROR: UndefVarError: `X2` not defined`.

Here, both X2 and Y2 becomes observables

julia> observed(osys)
2-element Vector{Equation}:
 X2(t) ~ 1 + X1(t)
 Y2(t) ~ 1 + Y1(t)

[ChrisRackauckas]

This test is passing on Catalyst? https://github.com/SciML/Catalyst.jl/actions/runs/9499687139/job/26181101288?pr=951#step:6:999

[TorkelE]

think so, yes