  We present a new computational framework for the Galerkin-collocation method
for double domain in the context of ADM 3+1 approach in numerical relativity.
This work enables us to perform high resolution calculations for initial sets
of two arbitrary black holes. We use the Bowen-York method for binary systems
and the puncture method to solve the Hamiltonian constraint. The nonlinear
numerical code solves the set of equations for the spectral modes using the
standard Newton-Raphson method, LU decomposition and Gaussian quadratures. We
show convergence of our code for the conformal factor and the ADM mass. Thus,
we display features of the conformal factor for different masses, spins and
linear momenta.
