Mapping of multiple lines such as the fine-structure emission from [C ii] (157.7 mu m), [O iii] (52 and 88.4 mu m), and rotational emission lines from CO are of particular interest for upcoming line intensity mapping (LIM) experiments at millimeter wavelengths, due to their brightness features. Several upcoming experiments aim to cover a broad range of scientific goals, from detecting signatures of the epoch of reionization to the physics of star formation and its role in galaxy evolution. In this paper, we develop a semianalytic approach to modeling line strengths as functions of the star formation rate (SFR) or infrared luminosity based on observations of local and high-z galaxies. This package, LIMpy (Line Intensity Mapping in Python), estimates the intensity and power spectra of [C ii], [O iii], and CO rotational transition lines up to the J levels (1-0) to (13-12) based both on analytic formalism and on simulations. We develop a relation among halo mass, SFR, and multiline intensities that permits us to construct a generic formula for the evolution of several line strengths up to z similar to 10. We implement a variety of star formation models and multiline luminosity relations to estimate the astrophysical uncertainties on the intensity power spectrum of these lines. As a demonstration, we predict the signal-to-noise ratio of [C ii] detection for an EoR-Spec-like instrument on the Fred Young Submillimeter Telescope. Furthermore, the ability to use any halo catalog allows the LIMpy code to be easily integrated into existing simulation pipelines, providing a flexible tool to study intensity mapping in the context of complex galaxy formation physics.
