We describe a Hartree ensemble method to approximately solve the Heisenberg equations for the phi (4) model in 1 + 1 dimensions. We compute the energies and number densities of the quantum particles described by the phi field and find that the particles initially thermalize with a Bose-Einstein distribution for the particle density. Gradually, however, the distribution changes towards classical equipartition. Using suitable initial conditions quantum thermalization is achieved much faster than the onset of this undesirable equipartition. We also show how the numerical efficiency of our method can be significantly improved.