We show that a serial array of N nonuniform, underdamped Josephson junctions coupled piezoelectrically to a nanoelectromechanical ͑NEM͒ oscillator results in phase locking of the junctions. Our approach is based on a semiclassical solution to a set of coupled differential equations that were generated by the Heisenberg operator equations, which in turn are based on a model Hamiltonian that includes the following effects: the charging and Josephson energies of the junctions, dissipation in the junctions, the effect of a dc bias current, an undamped simple harmonic oscillator ͑representing the NEM͒, and an interaction energy ͑due to the piezoelectric effect͒ between the NEM and the junctions. Phase locking of the junctions is signaled by a step in the current-voltage ͑I-V͒ curve. We find the phase-locked states are ͑neutrally͒ stable at the bottom and top of the step but not for bias currents in the middle of the step. Using harmonic balance, we are able to calculate an analytical expression for the location of the resonance step, v step , in the I-V curve. Because of the multistability of the underdamped junctions, it is possible, with a judicious choice of initial conditions and bias current, to set a desired number N a ഛ N of junctions on the resonance step, with N − N a junctions in the zero-voltage state. We are also able to show that, when N a junctions are in the phase-locked configuration, the time-averaged energy of the NEM oscillator scales like N a 2 .