We model the large scale late time universe as a Lambda-CDM cosmology driven by cosmological constant and perfect dust fluid. Our aim is to find new solutions in the matter and Lambda epoch consistent with inflationary initial conditions, namely that to the far past in the matter era the cosmology tends to a flat FLRW solution. We identify the moduli degrees of freedom that parametrize the flat Lambda-dust FLRW solution and then promote these moduli to slowly varying functions of the spatial coordinates and show how to solve the Einstein equations in a comoving gradient expansion, controlled by the cosmological constant length scale. Our initial conditions ensure that the approximation remains under control to the far past of the matter era, and to the far future of Lambda domination. The solution is fully non-perturbative in the amplitude of the metric deformation, and we explicitly construct it to fourth order in derivatives. A general Lambda-dust universe dominated by Lambda in the future is characterized by a 3-metric and a stress tensor (with positive trace) defined on the future conformal boundary. The new cosmologies with inflationary initial conditions are characterized only by the boundary 3-metric, the stress tensor being locally determined entirely in terms of that metric.

Assuming homogeneous isotropic Lambda-CDM cosmology allows Lambda, spatial curvature and dark matter density to be inferred from large scale structure observations such as supernovae. The purpose of this paper is to extend this to allow observations to measure or constrain inhomogeneity and anisotropy. We obtain the general inhomogeneous anisotropic Lambda-CDM solution which is locally asymptotic to an expanding de Sitter solution as a late time expansion using Starobinsky's method (analogous to the `holographic renormalization' technique in AdS/CFT) together with a resummation of the series. The dark matter is modeled as perfect dust fluid. The terms in the expansion systematically describe inhomogeneous and anisotropic deformations of an expanding FLRW solution, and are given as a spatial derivative expansion in terms of data characterizing the solution - a 3-metric and a perturbation of that 3-metric. Leading terms describe inhomogeneity and anisotropy on the scale set by the cosmological constant, approximately the horizon scale today. Higher terms in the expansion describe shorter scale variations. We compute the luminosity distance-redshift relation and argue that comparison with current and future observation would allow a partial reconstruction of the characterizing data. We also comment on smoothing these solutions noting that geometric flows (such as Ricci flow) applied to the characterizing data provide a canonical averaging method.

TeVeS uses a dynamical vector field with timelike unit norm constraint to specify a preferred local frame. When matter moves slowly in this frame - the so-called quasi-static regime - Modified Newtonian Dynamics (MoND) results. Theories with such vectors (such as Einstein-aether) are prone to the vector dynamics forming singularities which render their classical evolution problematic. Here we analyse the dynamics of the vector in TeVeS in various situations. We find that, quite generically, the vector field develops caustic singularities on time scales of order the gravitational in-fall time. Having shown singularity formation is generic with or without matter, Bekenstein's original formulation of TeVeS appears dynamically problematic. We argue that by modifying the vector field kinetic terms to the more general form used by Einstein-Aether this problem may be avoided.