In this paper, we describe the first steps towards fully non-perturbative cosmology. We explain why the conventional methods used by cosmologists based on the ADM formulation are generally inadequate for this purpose and why it is advantageous instead to adapt the harmonic formulation pioneered and utilized in mathematical and numerical relativity. Here we focus on using this approach to evaluating the linear mode stability in homogeneous and nearly homogeneous backgrounds and devising a valid scheme and diagnostics for numerical computation. We also briefly touch on the relevance of these methods for extracting cosmological observables from non-perturbative simulations.