A novel, compact mathematical formulation is presented to describe the dynamic rocking response of single and double block systems subjected to gravity and/or ground excitation. The derivation of the closed-form solutions for impact and motion is based on the Euler-Lagrange equation and the conservation of angular momentum, and combines all the different cases of possible block relative rotating and impact modes (16 in total) into a single set of equations without the need of transient expressions. The derived equations that describe the impact modes are the equivalent to the expression derived by Housner and depend on the angular velocity of the blocks before impact. The analytical model is integrated numerically via an ad hoc algorithm and its reliability & accuracy are verified after various self-consistency tests and comparisons with the literature. In addition, several shaking table experiments were conducted in EQUALS laboratory in Bristol, setup constructed to test free and forced rocking motion of single and double block configurations. The error margins of the measurements are determined, and the extracted data are in good agreement with the numerical results for most examined cases. The ideal Housner restitution coefficient of single block impact to a rigid base is adjusted to match experimental conditions, and it is found to be correlated with the block aspect ratio. The forced rocking of a two-block system is shown to exhibit numerous different response patterns depending on the excitation conditions. The integrated model is finally applied to produce normalised overturning maps for double block systems, subjected to singlepulse sine inputs, which uncover the existence of a fractal-type behaviour. This previously unsuspected trait of multi-block systems is reminiscent of the chaotic behaviour exhibited by a classical double pendulum and suggests that the risk of overturning can only be evaluated on a probabilistic sense.