Localization marks the breakdown of thermalization in subregions of quantum many-body systems in the presence of sufficiently large disorder. In this paper, we use numerical techniques to study thermalization and localization in a many-body system of coupled quantum harmonic oscillators obtained by discretizing a scalar quantum field theory in Minkowski spacetime. We consider a Gaussian initial state, constructed through a global mass quench, with a quadratic Hamiltonian, and solve for the system's exact dynamics without and with disorder in one and two spatial dimensions. We find that finite-size systems localize for sufficiently large disorder in both cases, such that the entanglement entropy of subregions retains its initial area-law behavior, and the system no longer develops long-range correlations. To probe the thermalization-to-localization transition further, we define a frequency gap ratio that measures adjacent gaps in the phase space eigenvalues of the Hamiltonian and study how it varies with disorder strength and system size. We find signatures of a chaotic regime at intermediate disorder in two spatial dimensions and argue that it is a finitesize effect, such that the system would localize for arbitrarily small disorder in the continuum in both one and two spatial dimensions, consistent with Anderson localization. Lastly, we use the frequency gap ratio to argue that in three spatial dimensions, on the other hand, the system would only localize for disorder strengths above a critical value in the continuum, again consistent with Anderson localization.