Non-linear evolution of the parametric instability of inertial waves inherent to eccentric discs is studied by way of a new local numerical model. Mode coupling of tidal deformation with the disc eccentricity is known to produce exponentially growing eccentricities at certain mean-motion resonances. However, the details of an efficient saturation mechanism balancing this growth still are not fully understood. This paper develops a local numerical model for an eccentric quasi-axisymmetric shearing box which generalizes the often-used Cartesian shearing box model. The numerical method is an overall second-order well-balanced finite volume method which maintains the stratified and oscillatory steady-state solution by construction. This implementation is employed to study the non-linear outcome of the parametric instability in eccentric discs with vertical structure. Stratification is found to constrain the perturbation energy near the mid-plane and localize the effective region of inertial wave breaking that sources turbulence. A saturated marginally sonic turbulent state results from the non-linear breaking of inertial waves and is subsequently unstable to large-scale axisymmetric zonal flow structures. This resulting limit-cycle behaviour reduces access to the eccentric energy source and prevents substantial transport of angular momentum radially through the disc. Still, the saturation of this parametric instability of inertial waves is shown to damp eccentricity on a time-scale of a thousand orbital periods. It may thus be a promising mechanism for intermittently regaining balance with the exponential growth of eccentricity from the eccentric Lindblad resonances and may also help explain the occurrence of ‘bursty’ dynamics such as the superhump phenomenon.