That paper needs an update. Taylor methods utilize AD effectively to get something that seems efficient. But there's a lot of theory to explain why they aren't efficient in practice, mostly because of truncation error coefficient sizes. Also, since the time of the article, there have been much better high order explicit RK methods that have been derived, and improved schemes for adaptive order explicit extrapolation. When we benchmark these against Taylor methods which use an optimal form of AD (provided by TaylorIntegration.jl), we still see about a 10x performance difference when at the edge of Float64 range:

https://benchmarks.juliadiffeq.org/html/DynamicalODE/Henon-H...

https://benchmarks.juliadiffeq.org/html/DynamicalODE/Quadrup...

As you get outside of Float64 range, the parallelized explicit extrapolation methods seem to do better. This is something we're going to be getting a lot more precise benchmarks in JuliaDiffEq now though, since we have optimized generic implementations of each of the methods and can thus use Float128 and arbitrary precision on all of the methods simultaneously to really figure out the cutoff.

But thinking about efficiency and accuracy misses the point of Taylor methods. The Taylor methods have a way of getting rigorous error bounds on the integration (like interval arithmetic, but for the continuous ODE), which is the real reason why they can be useful, with the extra cost of course.

High order means much less dependence on step size. Automatic (exact) differentiation means no roundoff from finite differences with small step size. They use a very clever and involved code generation approach, I have just hardcoded the recurrences in my implementation (much simpler). https://github.com/m4r35n357/ODE-Playground

It is not a black box approach like RK4, it depends on the model using a (hardly at all) restricted set of functions. Piecewise definitions are OK though.