Equatorial plasma bubbles (spread F/ionospheric irregularities) are plasma density irregularities or depletion in the ionosphere. It can be observed by radio waves, optical and in situ instruments during postsunset in equatorial and low‐latitude regions. Severe scintillations in radio waves can be caused by the presence of plasma bubbles. Therefore, studies of plasma bubbles become one of the most important and hot topics. Many studies have been devoted to the physical mechanism and characteristics of plasma bubbles in the ionosphere. However, some critical issues (e.g., day‐to‐day variability, complex structures, and ions composition) are still not well known. In this study, a new two‐dimensional model of multi‐ion (H+ , He+, N+ , O+ , N2+ , NO+ and O2+) plasma bubbles was developed. The photoionization, chemical reaction, and recombination processes of multiple ions have been added to the new model. The present model was used to simulate the nonlinear evolutions of multiple plasma bubbles in the equatorial region. Results show that the merging, disconnection and connection processes of plasma bubbles can be reproduced in this simulation. In addition, simulations first present morphological structures of molecular ions (NO+, O2+) of multiple plasma bubbles, which are significantly different from plasma bubbles. Molecular ion NO+ is manifested as plasma blobs. However, molecular ion O2+ is manifested as bubbles at low altitudes but blobs at high altitudes. The complex structures of multiple plasma bubbles and morphological structures of molecular ions in the present simulation are mostly consistent with observations.