In Julia, this is implemented with the number type of TaylorSeries.jl https://github.com/JuliaDiff/TaylorSeries.jl and captured in the integration package TaylorIntegration.jl https://github.com/PerezHz/TaylorIntegration.jl which is then available in the DifferentialEquations.jl common interface:

https://docs.sciml.ai/dev/modules/DiffEqDocs/solvers/ode_sol...

The implementation described here is missing a lot of the optimizations that TaylorIntegration.jl employs (it makes the method no longer simple), but after doing some of the extra stuff you can find that Taylor integration methods can be good if you need really high accuracy. Some examples of this are seen in some of the SciMLBenchmarks, particularly https://benchmarks.sciml.ai/html/DynamicalODE/Henon-Heiles_e... and https://benchmarks.sciml.ai/html/DynamicalODE/Quadrupole_bos... .

One of the nice things about Taylor integration methods is that they can be used to build interval solutions of ODEs which guarantee the solution. That's described in some detail in this video: https://www.youtube.com/watch?v=o1h7BUW04NI .