  Phase reduction is a commonly used techinque for analyzing stable
oscillators, particularly in studies concerning synchronization and phase lock
of a network of oscillators. In a widely used numerical approach for obtaining
phase reduction of a single oscillator, one needs to obtain the gradient of the
phase function, which essentially provides a linear approximation of isochrons.
In this paper, we extend the method for obtaining partial derivatives of the
phase function to arbitrary order, providing higher order approximations of
isochrons. In particular, our method in order 2 can be applied to the study of
dynamics of a stable oscillator subjected to stochastic perturbations, a topic
that will be discussed in a future paper. We use the Stuart-Landau oscillator
to illustrate the method in order 2.
