I am looking for a Julia package that can provide high-order Taylor expansion coefficients of the solution for a given ODE system and initial conditions using automatic differentiation. It does not need to solve the ODE. Ideally, the initial conditions should be input as intervals to allow for interval arithmetic.

(If you are trying to fit the solution to a polynomial, you might consider just evaluating the solution at a set of Chebyshev points and doing Chebyshev polynomial interpolation, for which multiple packages exist. This will have high accuracy over a whole interval, as opposed to Taylor series whose accuracy is concentrated around a single point.)

Iām not trying to solve the IVP; I only want to obtain the Taylor series coefficients using automatic differentiation. Each coefficient should be the Jacobian matrix of the previous coefficient multiplied by ( f ).

Do you want the taylor expansion w.r.t initial conditions or parameters? You can just use TaylorSeries.jl or TaylorDiff.jl with DifferentialEquations.jl for that. If you are looking for the Taylor expansion in time, the solution is already sol(t) a continuous higher order interpolating function, and `sol(t,Val{i})`

gives the ith derivative.