The behavior of many natural and engineered systems is determined by oscillatory phenomena for which the input-output relationship can be described using phase models. The use of such models significantly reduces the complexity of control design, and enables the application of powerful semi-analytical methods for optimal control synthesis. In this paper, we examine the optimal control of a collection of neuron oscillators described by phase models. In particular, we employ Pontryagin's maximum principle to formulate the optimal control problem as a boundary value problem, which we then solve using the homotopy perturbation method. This iterative optimization-free technique is promising for neural engineering applications that involve nonlinear oscillatory systems for which phase model representations are feasible.