On the example of the Poynting-Thomson-Zener rheological model for solids, which exhibits both dissipation and wave propagation -with nonlinear dispersion relation -, we introduce and investigate a finite difference numerical scheme. Our goal is to demonstrate its properties and to ease the computations in later applications for continuum thermodynamical problems. The key element is the positioning of the discretized quantities with shifts by half space and time steps with respect to each other. The arrangement is chosen according to the spacetime properties of the quantities and of the equations governing them. Numerical stability, dissipative error and dispersive error are analysed in detail. With the best settings found, the scheme is capable of making precise and fast predictions. Finally, the proposed scheme is compared to a commercial finite element software, COMSOL, which demonstrates essential differences even on the simplest -elastic -level of modelling.