Based on very recent and promising ideas, stemming from the use of bubbles, we discuss an algorithm for the numerical simulation of the cubic nonlinear Schrödinger equation with harmonic potential (cNLS) in any dimension, that could easily be extended to other polynomial nonlinearities. This algorithm consists in discretizing the initial function as a sum of modulated complex gaussian functions (the bubbles), each one having its own set of parameters, and then updating the parameters according to cNLS. Numerically, we solve exactly the linear part of the equation and use the Dirac-Frenkel-MacLachlan principle to approximate the nonlinear part. We then obtain a grid free algorithm in any dimension whose efficiency compared with spectral methods is illustrated by numerical examples.