  Mapping of multiple lines such as the fine-structure emission from [CII]
(157.7 $\mu \text{m}$), [OIII] (52 \& 88.4 $\mu \text{m}$), and rotational
emission lines from CO are of particular interest for upcoming line intensity
mapping (LIM) experiments at millimetre 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
semi-analytic approach to modelling line strengths as functions of the star
formation rate (SFR) or infrared (IR) luminosity based on observations of local
and high-z galaxies. This package, $\texttt{LIMpy}$ (Line Intensity Mapping in
Python), estimates the intensity and power spectra of [CII], [OIII], 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 multi-line intensities that permits us to construct a generic formula
for the evolution of several line strengths up to $z \sim 10$. We implement a
variety of star formation models and multi-line 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 [CII]
detection for an EoR-Spec-like instrument on the Fred Young Submillimeter
Telescope (FYST). Furthermore, the ability to use any halo catalogue allows the
$\texttt{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.
