  We have used the axisymmetric chemo-hydrodynamical code WALKIMYA-2D to
numerically model and reproduce the physical and CO emission properties of the
jet-driven outflow from the intermediate-mass protostar Cep E, which was
observed at $\sim 800$au resolution in the CO $J=2\to 1$ line with the IRAM
interferometer. Our simulations take into account the observational constraints
available on the physical structure of the protostellar envelope to provide
constraints on the dynamics of the inner protostellar environment from the
study of the outflow/jet propagation away from the launch region. WALKIMYA-2D
successfully reproduces the main qualitative and quantitative features of the
Cep E outflow and the jet kinematics, naturally accounting for their time
variability. Signatures of internal shocks are detected as knots along the jet.
In the early times of the ejection process, the young emitted knots interact
with the dense circumstellar envelope through high-velocity, dissociative
shocks, which strongly decrease the CO gas abundance in the jet. As time
proceeds, the knots propagate more smoothly through the envelope and
dissociative shocks disappear after $\sim 10^3$ yr. The distribution of CO
abundance along the jet shows that the latter bears memory of the early
dissociative phase in the course of its propagation. Analysis of the velocity
field shows that the jet material mainly consists of gas entrained from the
circumstellar envelope and accelerated away from the protostar at $700$ au
scale. As a result, the overall jet mass loss rate appears higher than the
actual mass ejection rate by a factor $\sim 3$. Numerical modeling of the Cep E
jet-driven outflow and comparison with the CO observations have allowed us to
peer into the outflow formation mechanism with unprecedented detail and to
retrieve the history of the mass-loss events that have shaped the outflow.
