Based on recent 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 in any dimension, which could be easily extended to other polynomial nonlinearities. For the linear part of the equation, the algorithm consists in discretizing the initial function as a sum of modulated complex functions, each one having its own set of parameters, and then updating the parameters exactly so that the modulated function remains a solution to the equation. When cubic interactions are introduced, the Dirac-Frenkel-MacLachlan principle is used to approximate the time evolution of parameters. We then obtain a grid-free algorithm in any dimension, and it is compared to a spectral method on numerical examples.