From the solution of a differential equation using Vern6, Vern7, Vern8 or Vern9, how can I extract the dense interpolant polynomial coefficients?

Of course I could evaluate the polynomial on a grid of size order+1 inside a step and then compute the interpolation polynomial using Lagrange’s method, but I do not want to evaluate the polynomial because the evaluation of that polynomial in Float64 arithmetics creates too much numerical noise for my application. I can solve this issue by either:

- using BigFloat numeric type (tested in Julia)
- evaluating the interpolation polynomial using compensated floating point arithmetics: https://www-pequan.lip6.fr/~jmc/polycopies/Compensation-horner.pdf (I tested this successfully in a C++ code with Verner’s RK877 dense interpolant computed in double precision)

I would like to try method 2 in Julia but for that I need access to the interpolant coefficients.

Thank you.