By employing the time-dependent variational principle and the versatile multi-D2 Davydov trial states, in combination with the Green’s function method, we study the dynamics of the Tavis–Cummings model and the Holstein–Tavis–Cummings model in the presence of diagonal disorder and cavity–qubit coupling disorder. For the Tavis–Cummings model, time evolution of the photon population, the optical absorption spectra, and the hetero-entanglement between the qubits and the cavity mode are calculated by using the Green’s function method to corroborate numerically exact results of Davydov’s Ansätze. For the Holstein–Tavis–Cummings model, only the latter is utilized to simulate effects of disorder on the photon population dynamics and the absorption spectra. We have demonstrated that the complementary employment of analytical and numerical methods permits uncovering a fairly comprehensive picture of a variety of complex behaviors in disordered multidimensional polaritonic cavity quantum electrodynamics systems.