•A new python code (ColRadPy) solves collisional radiative & ionization balance equations.•Application of ColRadPy to fusion relevant neutral molybdenum.•Many metastable states are often required to accurately model heavy species.•A new method to identify the most important metastable is described. The properties of emitting ions in a plasma provides both potential for plasma diagnostics and key information required for plasma modeling. Generalized collisional radiative theory provides a powerful tool for the modeling of low and moderately dense plasmas. A new Python program is presented that solves the collisional radiative and ionization balance equations for application to fusion, laboratory, and astrophysical plasmas. It produces generalized coefficients that can be easily imported into existing plasma modeling codes and spectral diagnostics. An overview of the code is presented, along with selected results for applications to high-Z plasma facing components.