A novel finite-difference time-domain formulation for the modeling of electromagnetic wave propagation in frequency-dispersive liquid crystals with random orientation is presented. The dispersion of the complex ordinary and extraordinary permittivities of nematic liquid crystals is described by a generalized model based on the sum of partial fractions. The proposed dispersive model encompasses traditional approaches, such as Drude, Drude–Lorentz, and modified-Lorentz functions; it can also capture arbitrary dispersion properties of experimentally characterized nematic mixtures via the vector fitting technique. The accuracy of the formulation is demonstrated in a series of benchmark examples in optical and terahertz frequencies.