Context. Linear polarization of spectral lines, commonly known as the Goldreich-Kylafis effect within the star formation community, is one of the most underutilized techniques for probing magnetic fields in the dense and cold interstellar medium. Aims. In this study, we implement linear polarization of molecular spectral lines into the multilevel, non-local thermodynamic equilibrium radiative transfer code PYRATE. Methods. Different modes of polarized radiation are treated individually, with separate optical depths computed for each polarization direction. Our implementation is valid in the so-called strong magnetic field limit and is exact for either a system satisfying the large-velocity-gradient approximation, and/or for any system with a uniform magnetic field. We benchmark our implementation against analytical results and provide tests for various limiting cases. Results. In agreement with previous theoretical results, we find that in the multilevel case the amount of fractional polarization decreases when compared to the two-level approximation, but this result is subject to the relative importance between radiative and collisional processes. Finally, we post-process an axially symmetric, nonideal magnetohydrodynamic chemo-dynamical simulation of a collapsing prestellar core and provide theoretical predictions regarding the shape (as a function of velocity) of the polarization fraction of CO during the early stages in the evolution of molecular clouds.