  We present a numerical method to compute quasiequilibrium configurations of
close binary neutron stars in the pre-coalescing stage. A hydrodynamical
treatment is performed under the assumption that the flow is either rigidly
rotating or irrotational. The latter state is technically more complicated to
treat than the former one (synchronized binary), but is expected to represent
fairly well the late evolutionary stages of a binary neutron star system. As
regards the gravitational field, an approximation of general relativity is
used, which amounts to solving five of the ten Einstein equations (conformally
flat spatial metric). The obtained system of partial differential equations is
solved by means of a multi-domain spectral method. Two spherical coordinate
systems are introduced, one centered on each star; this results in a precise
description of the stellar interiors. Thanks to the multi-domain approach, this
high precision is extended to the strong field regions. The computational
domain covers the whole space so that exact boundary conditions are set to
infinity. Extensive tests of the numerical code are performed, including
comparisons with recent analytical solutions. Finally a constant baryon number
sequence (evolutionary sequence) is presented in details for a polytropic
equation of state with gamma=2.
