The study of shrinking, expanding, and strongly interacting bubbles at high Reynolds numbers is of significant interest for micro- and nanotechnologies. One of such interests is related to self-propulsion of bubbles due to non-linear interaction of bubble shape modes. In the present study bubble dynamics in potential flow is considered. The boundary element method (BEM) which offers a low computational cost and provides an accurate representation of bubble surface is employed for studies. To accelerate computations and increase problem size the fast multipole method (FMM) and graphics processors (GPUs) are used. For mesh stabilization, which appears to be an issue, a new parametric spherical filter based on spherical harmonic expansion is developed and implemented. The dynamics of high order surface modes of bubble at free and forced bubble oscillations is studied.