A simple and efficient method for simulating dispersed bubble flow is provided. Instead of modeling the complex hydrodynamics of numerous small bubbles explicitly, the method approximates the average motion of these bubbles using a continuum multiphase solver. The subgrid interactions among bubbles are computed using a new stochastic solver. Using the proposed scheme, complex scenes with millions of bubbles can be simulated efficiently.