Two difficulties associated with the computations of thermal vibrational correlation functions are discussed. The first one is the lack of a well-behaved expression that is valid at both high-temperature and T → 0 K limits. Specifically, if the partition function and the propagator are considered separately, then thermal vibrational correlation functions may have an indeterminate form 0/0 in the limit T → 0 K. This difficulty is resolved when the partition function and the propagator are jointly considered in the harmonic approximation, which allows a problematic term that emanates from the zero-point energy to be canceled out, thereby producing a thermal correlation function with a determinate form in the T → 0 K limit. The second difficulty is related to the multivaluedness of the vibrational correlation function. We show numerically that an improper selection of branch leads to discontinuities in the computed correlation function and incorrect vibronic spectra. We propose a phase tracking procedure that ensures continuity of both real and imaginary parts of the correlation function to recover the correct spectra. We support our findings by simulating the UV-vis absorption spectra of pentacene at 4 K and benzene at 298 K. Both are found to be in good agreement with their experimental counterparts.