The star formation rate density of the Universe is known to peak around z ~ 1 – 2 and has declined to the present day. The fundamental physical processes shaping this evolution, however, remain uncertain. Constraints on the molecular gas mass density evolution are particularly relevant to understanding these processes because this gas is the immediate fuel for star formation. Most studies rely on blind CO searches in deep fields, which require hundreds of hours of telescope time and then use those detections to construct the CO luminosity function over the sampled cosmological volume. We instead use a new hybrid method by taking advantage of existing data from the Plateau de Bure High-z Blue Sequence Survey 2 (PHIBSS2) program. The observations targeted normal star-forming main sequence galaxies at the cosmic peak of star formation, selected from the 3D-HST/CANDELS fields. We do a systematic search in the data for serendipitous “secondary” sources and use these to measure the CO luminosity function for a variety of CO transitions and redshift ranges. With these data, we derive the molecular gas mass density evolution over the redshift range z ~ 0.3 – 4, and find that our result is consistent with previous observations. We also add constraints on CO luminosity functions for previously unexplored transitions and redshift ranges such as CO(3– 2) at z ~ 1.5, CO(4–3) at z ~ 2.2, and CO(5–4) at z ~ 3.3. Our method can be a powerful new way of measuring the CO luminosity function and molecular gas mass density over cosmic time by leveraging existing data as well as broader bandwidth and higher sensitivity data from future facilities such as the ngVLA.