We study the time correlation functions of coupled linear Langevin dynamics with and without inertial effects, both analytically and numerically. The model equation represents the physical behavior of a harmonic oscillator in two or three dimensions in the presence of friction, additive noise, and an external field with both rotational and deformational components. This simple model plays pivotal roles in understanding more complicated processes. The analytical solution presented serves as a test of numerical integration schemes, and its derivation is presented in a fashion that allows it to be repeated directly in a classroom. While the results in the absence of fields (equilibrium) or confinement (free particle) are omnipresent in the literature, we write down, apparently for the first time, the full nonequilibrium results that may correspond, e.g., to a Hookean dumbbell embedded in a macroscopically homogeneous shear or mixed flow field. We demonstrate how the inertial results reduce to their noninertial counterparts in the nontrivial limit of vanishing mass. While the results are derived using basic integrations over Dirac delta distributions, we also provide alternative approaches involving (i) Fourier transforms, which seem advantageous only if the measured quantities also reside in Fourier space, and (ii) a Fokker--Planck equation and the moments of the probability distribution. The results, verified by numerical experiments, provide additional means of measuring the performance of numerical methods for such systems. It should be emphasized that this article provides specific details regarding the derivations of the time correlation functions as well as the implementations of various numerical methods, so that it can serve as a standalone piece for lessons in the framework of Itô stochastic differential equations and calculus.