  We present a numerically tractable method to solve exactly the evolution of a
N boson system with binary interactions. The density operator of the system rho
is obtained as the stochastic average of particular operators |Psi_1><Psi_2| of
the system. The states |Psi_{1,2}> are either Fock states |N:phi_{1,2}> or
coherent states |coh:phi_{1,2}> with each particle in the state phi_{1,2}. We
determine the conditions on the evolution of phi_{1,2} -which involves a
stochastic element- under which we recover the exact evolution of rho. We
discuss various possible implementations of these conditions. The well known
positive P-representation arises as a particular case of the coherent state
ansatz. We treat numerically two examples: a two-mode system and a
one-dimensional harmonically confined gas. These examples, together with an
analytical estimate of the noise, show that the Fock state ansatz is the most
promising one in terms of precision and stability of the numerical solution.
