  A new numerical method to construct binary black hole/neutron star initial
data is presented. The method uses three spherical coordinate patches; Two of
these are centered at the binary compact objects and cover a neighborhood of
each object; the third patch extends to the asymptotic region. As in the
Komatsu-Eriguchi-Hachisu method, nonlinear elliptic field equations are
decomposed into a flat space Laplacian and a remaining nonlinear expression
that serves in each iteration as an effective source. The equations are solved
iteratively, integrating a Green's function against the effective source at
each iteration. Detailed convergence tests for the essential part of the code
are performed for a few types of selected Green's functions to treat different
boundary conditions. Numerical computation of the gravitational potential of a
fluid source, and a toy model for a binary black hole field are carefully
calibrated with the analytic solutions to examine accuracy and convergence of
the new code. As an example of the application of the code, an initial data set
for binary black holes in the Isenberg-Wilson-Mathews formulation is presented,
in which the apparent horizons are located using a method described in Appendix
A.
