Abstract. The six Painlevé transcendents P I -P V I have both applications and analytic properties that make them stand out from most other classes of special functions. Although they have been the subject of extensive theoretical investigations for about a century, they still have a reputation for being numerically challenging. In particular, their extensive pole fields in the complex plane have often been perceived as 'numerical mine fields'. In the present work, we note that the Painlevé property in fact provides the opportunity for very fast and accurate numerical solutions throughout such fields. When combining a Taylor/Padé-based ODE initial value solver for the pole fields with a boundary value solver for smooth regions, numerical solutions become available across the full complex plane. We focus here on the numerical methodology, and illustrate it for the P I equation. In later studies, we will concentrate on mathematical aspects of both the P I and the higher Painlevé transcendents.Key words. Painlevé transcendents, P I equation, Taylor series method, Padé approximation, Chebyshev collocation method.1. Introduction. The six Painlevé transcendents were discovered about 100 years ago [4], [13], [23], [24], and have since received extensive attention. They satisfy nonlinear second order ODEs, cannot be reduced to other types special functions, and obey the Painlevé property of all its solutions being free from movable branch points (however, the solutions may have movable poles or movable isolated essential singularities) [18], [21]. They arise in many branches of physics, including classical and quantum spin models, scattering theory, self-similar solutions to nonlinear dispersive wave PDEs, string theory and random matrix theory. For a summary of their applications, with extensive references, we refer to Sections 32.13 -32.16 in the NIST Handbook of Mathematical Functions [21].The attention that the Painlevé functions have received over the last century has primarily been limited to their analytic properties and applications. The numerical computation of these functions, by contrast, has received comparatively little attention. With the exceptions of [8], [20], [22] (discussed later), numerical studies have so far mostly focused on pole-free intervals along the real line [17], [21], [25] or on a polefree region in the complex plane [11]. A possible explanation for this state of affairs is that the large pole fields in the complex plane associated with the Painlevé ODEs may be thought of as a challenge to standard numerical methods for the integration of ODEs.We argue in the present work that these pole fields provide, in fact, especially advantageous integration opportunities. A 'pole friendly' algorithm based on Padé approximation is devised that is capable of integrating accurately and stably through pole fields for distances in the 10 5 − 10 6 range. This allows numerical studies not only of 'near-fields', but also of transition processes towards 'far fields', for which an extensive body of asymptotic...