A numerical method is presented for simulating charged colloidal dispersions in electrolyte solutions. Utilizing a smoothed profile for colloid-solvent boundaries, efficient mesoscopic simulations are enabled for modeling dispersions of many colloidal particles exhibiting many-body electrostatic interactions. The validity of the method was examined for simple colloid geometries, and the efficiency was demonstrated by calculating stable structures of two-dimensional dispersions, which resulted in the formation of colloidal crystals.