Abstract. A Taylor method for solving an ordinary differential equation initial-value problem [Formula: see text], [Formula: see text], computes the Taylor series (TS) of the solution at the current point, truncated to some order, and then advances to the next point by summing the TS with a suitable stepsize. A standard ODE method (e.g., Runge–Kutta) treats function [Formula: see text] as a black box, but a Taylor solver requires [Formula: see text] to be preprocessed into a code-list of elementary operations that it interprets as operations on (truncated) TS. The trade-off for this extra work includes arbitrary order, typically enabling much larger stepsizes. For a standard function, such as [Formula: see text], this means evaluating [Formula: see text], where [Formula: see text] are TS. The sub-ODE method applies the ODE [Formula: see text], obeyed by [Formula: see text], to in-line this operation as [Formula: see text]. This gives economy of implementation: each function that satisfies a simple ODE goes into the “Taylor library” with a few lines of code—not needing a separate recurrence relation, which is the typical approach. Mathematically, however, the use of sub-ODEs generally transforms the original ODE into a differential-algebraic system, making it nontrivial to ensure a sound system of recurrences for Taylor coefficients. We prove that, regardless of how many sub-ODEs are incorporated into [Formula: see text], this approach guarantees a sound system. We introduce our sub-ODE-based MATLAB ODE solver and show that its performance compares favorably with solvers from the MATLAB ODE suite.