  Modelling the molecular gas that is routinely detected through CO
observations of high-redshift galaxies constitutes a major challenge for ab
initio simulations of galaxy formation. We carry out a suite of cosmological
hydrodynamic simulations to compare three approximate methods that have been
used in the literature to track the formation and evolution of the simplest and
most abundant molecule, H$_2$. Namely, we consider: i) a semi-empirical
procedure that associates H$_2$ to dark-matter haloes based on a series of
scaling relations inferred from observations, ii) a model that assumes chemical
equilibrium between the H$_2$ formation and destruction rates, and iii) a model
that fully solves the out-of-equilibrium rate equations and accounts for the
unresolved structure of molecular clouds. We study the impact of finite spatial
resolution and show that robust H$_2$ masses at redshift $z\approx 4$ can only
be obtained for galaxies that are sufficiently metal enriched in which H$_2$
formation is fast. This corresponds to H$_2$ reservoirs with masses
$M_{\mathrm{H_2}}\gtrsim 6\times 10^9 \mathrm{M}_\odot$. In this range,
equilibrium and non-equilibrium models predict similar molecular masses (but
different galaxy morphologies) while the semi-empirical method produces less
H$_2$. The star formation rates as well as the stellar and H$_2$ masses of the
simulated galaxies are in line with those observed in actual galaxies at
similar redshifts that are not massive starbursts. The H$_2$ mass functions
extracted from the simulations at $z\approx 4$ agree well with recent
observations that only sample the high-mass end. However, our results indicate
that most molecular material at high $z$ lies yet undetected in reservoirs with
$10^9<M_{\mathrm H_2}<10^{10} \mathrm{M}_\odot$.
