At a time when galaxy surveys and other observations are reaching unprecedented sky coverage and precision, it seems timely to investigate the effects of general relativistic nonlinear dynamics on the growth of structures and on observations. Analytic inhomogeneous cosmological models are an indispensable way of investigating and understanding these effects in a simplified context. In this paper, we develop exact inhomogeneous solutions of general relativity with pressureless matter (dust, describing cold dark matter) and cosmological constant Λ, which can be used to model an arbitrary initial matter distribution along one line of sight. In particular, we consider the second class Szekeres models with Λ, and split their dynamics into a flat ΛCDM background and exact nonlinear inhomogeneities, obtaining several new results. One single metric function Z describes the deviation from the background. We show that F, the time dependent part of Z, satisfies the familiar linear differential equation for δ, the first-order density perturbation of dust, with the usual growing and decaying modes. In the limit of small perturbations, δ≈F as expected, and the growth of inhomogeneities links up exactly with standard perturbation theory. In particular, we exhibit an exact conserved curvature variable, necessary for the existence of the growing mode, which is the nonlinear extension of the first-order curvature perturbation. We provide analytic expressions for the exact nonlinear δ and the growth factor in our models. For the case of over-densities, we find that, depending on the initial conditions, the growing mode may or may not lead to a pancake singularity, analogous to a Zel’dovich pancake. This is in contrast with the Λ=0 pure Einstein-de Sitter background where, at any given point in comoving (Lagrangian) coordinates, pancakes will always occur. Analyzing the covariant variables associated with the space-time, we derive the associated dynamical system, which we are able to decouple and reduce to two differential equations, one for ΩΛ representing the background dynamics, and one for δ describing the dynamics of the inhomogeneities. Our models are Petrov type D, which we show by explicitly deriving the only nonzero Weyl scalar Ψ2, which does not depend on Λ. Since this is the only Weyl contribution to the geodesic deviation equation, Λ can only contribute to lensing through its contribution to the background expansion.