FBDF provides bad interpolant for `dy/dt` in `dy ~ sin(t)`

[topolarity]

Describe the bug 🐞

FBDF gives a very bad interpolant for sol(t, Val{1}) of a simple ODE system dy ~ sin(t)

Minimal Reproducible Example 👇

julia> using OrdinaryDiffEq, Plots, LinearAlgebra
julia> ode_f1(u, _, t) = sin(t)
julia> ode_f2(u, _, t) = [u[2], u[2] - sin(t)]
julia> mass_matrix = Diagonal([1.0,0.0])
julia> prob1 = ODEProblem(ode_f1, 0.0, (0., 2π), nothing, abstol=1e-6, reltol=5e-2);
julia> prob2 = ODEProblem(ODEFunction(ode_f2; mass_matrix), [0.0, 0.0], (0., 2π), nothing, abstol=1e-6, reltol=5e-2);
julia> sol1 = solve(prob1, FBDF());
julia> sol2 = solve(prob2, FBDF());

julia> t′ = 0.0:1e-3:2π
julia> plot(t′, [u[1] for u in sol1(t′, Val{1}; idxs=1)]; title="FBDF (ODE)", label="sol1(t, Val{1})")
julia> plot(t′, hcat([[sol2(t, Val{1}; idxs=1), sol2(t, Val{0}; idxs=2)] for t in t′]...)'; label=["sol2(t, Val{1})" "sol2(t, Val{0}) (dy = x)"], title="FBDF (DAE)")

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

(tmp) pkg> status
Status `~/repos/tmp/tmp/Project.toml`
  [459566f4] DiffEqCallbacks v2.36.1
  [1dea7af3] OrdinaryDiffEq v6.70.0
  [c3572dad] Sundials v4.23.2

[ChrisRackauckas]

It's essentially the same problem as https://github.com/SciML/Sundials.jl/issues/452, though in this case we do actually have the possibility to define the corrected interpolation function. This is what I discussed in one of the chats as a good next direction with the BDFs.