Rodas5(P) 2nd and third derivatives

[ChrisRackauckas]

Fixes #2101

Second derivative looks fine, third derivative looks weird.

using OrdinaryDiffEq, ODEProblemLibrary
prob = prob_ode_linear
sol1 = solve(prob, Rodas5P(), dt = 1 // 2^(8), dense = true)
sol2 = solve(prob, Vern9(), dt = 1 // 2^(8), dense = true)

using Plots
t = 0.0:0.01:1.0
plot(t, Array(sol1(t, Val{1})))
plot!(t, Array(sol2(t, Val{1})))

plot(t, Array(sol1(t, Val{2})))
plot!(t, Array(sol2(t, Val{2})))

plot(t, Array(sol1(t, Val{3})))
plot!(t, Array(sol2(t, Val{3})))

[gstein3m]

Here is another example which could explain the behaviour:

using OrdinaryDiffEq, Plots

const n = 5

function f(du,u,p,t)
    du[1] = n*t^(n-1)
end
tspan = [0.0, 1.0]; u0 = [0.0]
prob1 = ODEProblem(ODEFunction(f), u0, tspan)
  
sol1 = solve(prob1, Rodas5(), dt = 1 // 2^(8), dense = true)
println(sol1.destats)
err_max = maximum(abs.(sol1[1,:] .- sol1.t.^n))
println("err_max:",err_max)
sol2 = solve(prob1, Vern9(), dt = 1 // 2^(8), dense = true)

t = 0.0:0.01:1.0
plot(sol1(t, Val{3}))
plot!(sol2(t, Val{3}))

The polynimial of order 5 can be integrated exactly by Rodas5 / 5P.
Since the interpolation error is not checked, only 5 time steps are taken regardless of the choice of dt. The thrid derivative of the fourth order interpolation is a linear function. This cannot give a good approximation of a second order polynomial within only 5 steps.
I therefore think that the proposed implementation is correct.

Independently of this, the evaluation of the first derivative for Rodas4/4P/4P2 does not work for me.
However, this may also be due to my suggested changes in PR #2092

[ChrisRackauckas]

I see. We should update that test case sometime then, but for now it does seem correct.