Ionic Bose polarons are quantum entities emerging from the interaction between an ion and a Bose-Einstein condensate (BEC), featuring long-ranged interactions that can compete with the gas healing length. This can result in strong interparticle correlations and enhancement of gas density around the ion. One possible approach to describe this complex system with high accuracy relies on numerical treatment such as the quantum Monte Carlo (QMC) techniques. Nevertheless, it is computationally very expensive and does not easily allow to study the system dynamics. On the other hand, a mean-field based variational ansatz in the co-moving frame can capture a sizeable change in the gas density. We apply it to the case of regularized ion-atom potential and find that it qualitatively reproduces the full numerical results. In addition, we also study the system of two pinned ions, focusing on their effective interaction induced by the bath. This approach seems to be promising for studying transport and nonequilibrium dynamics of charged (bi)polarons in condensed media.