We describe an extension of the Taylor method for the numerical solution of ODEs that uses Padé approximants to obtain extremely precise numerical results. The accuracy of the results is essentially limited only by the computer time and memory, provided that one works in arbitrary precision. In this method the stepsize is adjusted to achieve the desired accuracy (variable stepsize), while the order of the Taylor expansion can be either fixed or changed at each iteration (variable order). As an application, we have calculated the periodic solutions (limit cycle) of the van der Pol equation with an unprecedented accuracy for a large set of couplings (well beyond the values currently found in the literature) and we have used these numerical results to validate the asymptotic behavior of the period, of the amplitude and of the Lyapunov exponent reported in the literature. We have also used the numerical results to infer the formulas for the asymptotic behavior of the fast component of the period and of the maximum velocity, which have never been calculated before. •We introduce a Pade–Taylor method suitable for stiff ODEs.•Can obtain very accurate numerical solutions.•We have tested the method on the van der Pol equation (stiff and non-stiff).•We have obtained the most precise estimates of period, amplitude, Lyapunov exponents, etc.•We have used Richardson and Richardson–Pade extrapolation to derive the asymptotic behaviors.