Rocks exhibit beyond-Hookean, delayed and damped elastic, behaviour (creep, relaxation etc.). In many cases, the Poynting–Thomson–Zener (PTZ) rheological model proves to describe these phenomena successfully. A forecast of the PTZ model is that the dynamic elasticity coefficients are larger than the static (slow-limit) counterparts. This prediction has recently been confirmed on a large variety of rock types. Correspondingly, according to the model, the speed of wave propagation depends on frequency, the high-frequency limit being larger than the low-frequency limit. This frequency dependence can have a considerable influence on the evaluation of various wave-based measurement methods of rock mechanics. As experience shows, commercial finite element softwares are not able to properly describe wave propagation, even for the Hooke model and simple specimen geometries, the seminal numerical artefacts being instability, dissipation error and dispersion error, respectively. This has motivated research on developing reliable numerical methods, which amalgamate beneficial properties of symplectic schemes, their thermodynamically consistent generalization (including contact geometry), and spacetime aspects. The present work reports on new results obtained by such a numerical scheme, on wave propagation according to the PTZ model, in one space dimension. The simulation outcomes coincide nicely with the theoretically obtained phase velocity prediction.