We present a multi-dimensional continuum mathematical model for modeling the growth of a symbiotic biofilm system. We take a dual-species namely, the Streptococcus–Veillonella sp. biofilm system as an example for numerical investigations. The presented model describes both the cooperation and competition between these species of bacteria. The coupled partial differential equations are solved by using an integrative finite element numerical strategy. Numerical examples are carried out for studying the evolution and distribution of the bio-components. The results demonstrate that the presented model is capable of describing the symbiotic behavior of the biofilm system. However, homogenized numerical solutions are observed locally. To study the homogenization behavior of the model, numerical investigations regarding on how random initial biomass distribution influences the homogenization process are carried out. We found that a smaller correlation length of the initial biomass distribution leads to faster homogenization of the solution globally, however, shows more fluctuated biomass profiles along the biofilm thickness direction. More realistic scenarios with bacteria in patches are also investigated numerically in this study.