Symplectic numerical schemes for reversible dynamical systems predict the solution reliably over large times as well, and are a good starting point for extension to schemes for simulating irreversible situations like viscoelastic wave propagation and heat conduction coupled via thermal expansion occuring in rocks, plastics, biological samples etc. Dissipation error (artificial nonpreservation of energies and amplitudes) of the numerical solution should be as small as possible since it should not be confused with the real dissipation occurring in the irreversible system. In addition, the other well-known numerical artefact, dispersion error (artificial oscillations emerging at sharp changes), should also be minimal to avoid confusion with the true wavy behavior. The continuum thermodynamical aspects (respect for balances with fluxes, systematic constitutive relationships between intensive quantities and fluxes, the second law of thermodynamics with positive definite entropy production, and the spacetime-based kinematic viewpoint) prove valuable for obtaining such extended schemes and for monitoring the solutions. Generalizing earlier works in this direction, here, we establish and investigate such a numerical scheme for one-dimensional viscoelastic wave propagation in the presence of heat conduction coupled via thermal expansion, demonstrating long-term reliability and the applicability of thermodynamics-based quantities in supervising the quality of the solution.