After an historical review of several approaches to predict the spatial and energy distributions of ions confined in a Paul trap and subjected to collisions with neutrals, the paper focuses on the temporal invariance method. Two improvements of this method are studied. The first one is used to probe the shape of the equilibrium distribution around a given value of the dynamic state of the ion. The second one is an enhancement to compute the solution by the Monte Carlo technique: a large increase of the performances of the variance reduction is achieved for light buffer gases (He and Ne) by using a control model of the tested device.The application to the computation of the equilibrium distribution of Cs+ ions subjected to collisions with He, Ne and Ar shows that the distribution is not Gaussian, even with the lightest buffer gas atom, and the magnitude of the tail of the distribution dramatically increases with the collisional atom mass. In addition, the possibility to extend the studies to other confinement devices is discussed if the efficiency obtained by means of the use of a control model is increased.