Silent and wrong derivatives of observed variables

[hersle]

Derivatives of observed variables calculated with ODESolution(t, Val{1}, ...) are incorrect. On v9.12.2 (after #2574), the simple test

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

@variables x(t) y(t)
@named sys = ODESystem([D(x) ~ 1, y ~ x^2], t)
sys = structural_simplify(sys) # move y to observed variables
prob = ODEProblem(sys, [x => 0.0], (0.0, 1.0)) # analytical solution: x = t, y = t^2
sol = solve(prob)

y′_obs(t) = sol(t, Val{1}, idxs=y) # request derivative of observed variable y
y′_anal(t) = 2 * t # analytical solution y′ = 2*x*x′ = 2*t
@test y′_obs(1.0) ≈ y′_anal(1.0)

fails with

Test Failed at bug.jl:14
  Expression: y′_obs(1.0) ≈ y′_anal(1.0)
   Evaluated: 0.9999999999999751 ≈ 2.0

Is this the general case mentioned in #2574? Having this work generally would be both very convenient and a natural extension to the solution calling interface!

Also, the current behavior is to return the wrong derivative with no error or warning. I think this is severe and likely to cause hard-to-find bugs in user code. Is it possible to identify the cases when y is an observed variable and the methods of #2574 do not apply? In those cases, for example, could ODESolution(t, Val{N >= 1}, idxs=y) issue an error or warning?