Phase resetting is a common experimental approach to investigating the behaviour of oscillating neurons. Assuming repeated spiking or bursting, a phase reset amounts to a brief perturbation that causes a shift in the phase of this periodic motion. The observed effects not only depend on the strength of the perturbation, but also on the phase at which it is applied. The relationship between the change in phase after the perturbation and the unperturbed old phase, the so-called phase resetting curve, provides information about the type of neuronal behaviour, although not all effects of the nature of the perturbation are well understood. In this chapter, we present a numerical method based on the continuation of a multi-segment boundary value problem that computes phase resetting curves in ODE models. Our method is able to deal effectively with phase sensitivity of a system, meaning that it is able to handle extreme variations in the phase resetting curve, including resets that are seemingly discontinuous. We illustrate the algorithm with two examples of planar systems, where we also demonstrate how qualitative changes of a phase resetting curve can be characterised and understood. A seven-dimensional example emphasises that our method is not restricted to planar systems, and illustrates how we can also deal with non-instantaneous, time-varying perturbations.