The complex behavior of heterogeneous cell populations can be accurately described by a class of mathematical models widely known as cell population balance models. They are first order partial integro-differential equations, whose formulation is fully defined by the single-cell reaction and division rates and the partition probability density function for the description of the mechanism which defines the distribution of cellular material amongst the daughter cells at cell division. Numerical algorithms have been developed for the numerical solution of cell population balance models due to the difficulty in deriving analytical solutions. However, their efficient solution still remains a challenging task; the main challenge comes from the fact that the boundaries of the intracellular state space are typically not known a priori. A novel free boundary algorithm has been recently developed treating with success the challenges of unknown boundaries and the co-existence of solutions differing by orders of magnitude, based on a transformation of the independent variable with respect to the mean expression level of the intracellular content. We present an extension of this free boundary algorithm in systems with two intracellular species. As a test bed of our algorithm we considered the case of linear single-cell reaction rates with respect to both species. Equal partitioning of cellular material at cell division exhibits oscillatory behavior for the number density function normalized around the mean intracellular content. The period of oscillations depends on the sharpness of the division mechanism, while unequal partitioning at cell division destroys the oscillations. We also present a parametric study of steady state solutions in E-Coli populations with a genetic toggle switch. Our target is to elucidate the effect of heterogeneity in the extent of the bistability region which appears at the single-cell level.