  We propose an efficient stochastic method to implement numerically the
Bogolubov approach to study finite-temperature Bose-Einstein condensates. Our
method is based on the Wigner representation of the density matrix describing
the non condensed modes and a Brownian motion simulation to sample the Wigner
distribution at thermal equilibrium. Allowing to sample any density operator
Gaussian in the field variables, our method is very general and it applies both
to the Bogolubov and to the Hartree-Fock Bogolubov approach, in the equilibrium
case as well as in the time-dependent case. We think that our method can be
useful to study trapped Bose-Einstein condensates in two or three spatial
dimensions without rotational symmetry properties, as in the case of
condensates with vortices, where the traditional Bogolubov approach is
difficult to implement numerically due to the need to diagonalize very big
matrices.
