We present a new method for simulating bubble dynamics called Level Set Indirect Boundary Element Method (LSBEM). This method tries to combine the advantages of LSM with BEM. As we know, the Level Set Method (LSM), which is a long-familiar computational technique used for tracking a propagating interface over time, has the strength in accurately handling topological complexities and changes. The Boundary Element Method (BEM), on the other hand, is known to serve well in conserving computational effort by reducing the dimensions of the problem by one. The novelty of this work is that while keeping this advantage of BEM, LSBEM simplifies the representation of the interface of multi-bubbles by using LSM. Thus advantages from both methods are conserved. A number of techniques are applied to ensure solution convergence and numerical accuracy. For instance, effort is made to avoid singularities in calculation by defining two sets of source and control points on the mesh which will never overlap one another; also work is done to ensure solution accuracy by reinitializing the level set function using Fast Marching Method (FMM) after every two timesteps.