We report the results of three-dimensional Newtonian calculations of neutron-star binary coalescence using smooth particle hydrodynamics (SPH). Using a relaxation technique, we construct hydrostatic equilibrium models of close neutron-star binaries in synchronized circular orbits. We use a simple polytropic equation of state with GAMMA = 2 to represent cold nuclear matter, and we assume that the mass ratio q = 1, as observed in all known neutron-star binary systems. Using SPH, we study the dynamical stability of these hydrostatic equilibrium models. In a sequence of models with decreasing binary separation we find that dynamical instability sets in slightly before the point along the sequence where the surfaces of the two stars come into contact. This is in agreement with the known stability properties of the solutions of the classical Darwin problem for two identical, incompressible components. We find that the initial stage of the instability, consisting in the steady merging of the two stars into a single ellipsoidal object, is completed in about one orbital period. At this point sudden mass shedding is triggered, resulting in the rapid removal of matter from the central object through two outgoing spiral arms. This results in the rapid redistribution of matter in the system until a new, nearly axisymmetric, differentially rotating equilibrium structure has formed. Using the quadrupole approximation, we follow the emission of gravitational radiation from the onset of dynamical instability to the establishment of axial symmetry. To support our results, we present several test-bed calculations which use SPH for binary systems. We consider axisymmetric, head-on collisions between two identical GAMMA = 2 polytropes and compare our SPH results to those of previous finite-difference calculations. Most importantly, we calculate solutions of the Roche and Darwin problems for polytropes with a wide range of adiabatic indices, 5/3 less-than-or-equal-to GAMMA less-than-or-equal-to 10. We find good agreement with known analytical results, in both the nearly incompressible and highly compressible limiting regimes. These calculations provide stringent tests of our method's ability to hold stable binaries in equi-librium and to identify terminal points or the onset of dynamical instability along equilibrium sequences of close binaries. Such tests are crucial for establishing the credibility of numerical results and, in particular, of computed gravitational radiation waveforms.