  Pulsar timing arrays (PTAs) might detect gravitational waves (GWs) from
massive black hole (MBH) binaries within this decade. The signal is expected to
be an incoherent superposition of several nearly-monochromatic waves of
different strength. The brightest sources might be individually resolved, and
the overall deconvolved, at least partially, in its individual components. In
this paper we extend the maximum-likelihood based method developed in Babak &
Sesana 2012, to search for individual MBH binaries in PTA 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 multi-search genetic algorithm to maximize the likelihood function
and look for sources in synthetic datasets. On datasets characterized by white
Gaussian noise plus few injected sources with signal-to-noise ratio (SNR) in
the range 10-60, our search algorithm performs well, recovering all the
injections with no false positives. Individual source SNRs are estimated within
few % of the injected values, sky locations are recovered within few degrees,
and frequencies are determined with sub-Fourier bin precision.
