Pulsar timing arrays might detect gravitational waves from massive black hole binaries within this decade. The signal is expected to be an incoherent superposition of several nearly monochromatic waves of different strengths. The brightest sources might be individually resolved, and the overall deconvolved, at least partially, in their individual components. In this paper we extend the maximum-likelihood-based method developed in Babak and Sesana [Phys. Rev. D 85, 044034 (2012)] to search for individual massive black hole binaries in pulsar timing array data. We model the signal as a collection of circular monochromatic binaries, each characterized by three free parameters: two angles defining the sky location and the frequency. We marginalize over all other source parameters, and we apply an efficient multisearch genetic algorithm to maximize the likelihood function and look for sources in synthetic data sets. On data sets characterized by white Gaussian noise plus few injected sources with signal-to-noise ratio in the range 10--60, our search algorithm performs well, recovering all the injections with no false positives. Individual source signal-to-noise ratios are estimated within a few percent of the injected values, sky locations are recovered within a few degrees, and frequencies are determined with sub-Fourier-bin precision.