We obtain an elegant and useful description of the dynamics of Szekeres dust models (in their full generality) by means of `quasi-local' scalar variables constructed by suitable integral distributions that can be interpreted as weighed proper volume averages of the local covariant scalars. In terms of these variables, the field equations and basic physical and geometric quantities are formally identical to their corresponding expressions in the spherically symmetric LTB dust models. Since we can map every Szekeres model to a unique LTB model, rigorous results valid for the latter models can be readily generalized to a non-spherical Szekeres geometry. The new variables lead naturally to an initial value formulation in which all scalars are expressed as scaling laws in terms of their values at an arbitrary initial space slice. These variables also yield a significant simplification of numerical work, since the fluid flow evolution equations become a system of autonomous ordinary differential equations subjected to algebraic constraints containing the information on the deviations from spherical symmetry. As an example of how this formalism can be applied, we show that spherical symmetry is stable against small dipole-like perturbations. This new approach to the dynamics of the Szekeres solutions has an enormous potential for dealing with a wide variety of theoretical issues and for constructing non-spherical models of cosmological inhomogeneities to fit observational data.