Abstract <p>The paper sets forth mathematical modeling, a method, and an algorithm for numerically solving the problem of wave heat transfer in a nonequilibrium nonlinear bounded medium. The modeling of wave heat transfer is based on the nonlinear Vernotte–Cattaneo–Lykov law, which is affected by isolated rectangular temperature pulses periodically applied to the boundary. The nonlinearity is due to the temperature dependence of the thermal conductivity, which is modeled as a power function with a non-integer exponent. This results in a nonlinear system of difference equations for the nodal temperature values on the upper time layer (implicit approximation). To avoid using iterative methods for solving this system and to employ an efficient Thomas algorithm method, the temperature dependencies of thermal conductivity on the upper time layer are linearized at each time layer using linear time extrapolation from the two previous time layers (where the temperature distribution is already known) to the upper time layer. The algorithm that we developed produced the results of numerical experiments, which allowed us to analyze the dynamic and kinematic characteristics of thermal waves, the speeds of movement and reflection from the opposite wall of fronts in the form of finite discontinuities in the temperature distribution. New phenomena were discovered, including the formation of thermal shock waves due to the superposition on previous pulses of subsequent thermal pulses moving through the heated space at a higher speed. Wave heat transfer in a nonlinear space was found to be similar to heat transfer in a linear space provided it had sources of thermal energy.</p>