A numerical procedure for the calculation of space distribution of the thermal field in thermally inhomogeneous solids with thermal memory induced by a harmonically modulated surface heat source is presented. The procedure is based on the division of the complex hyperbolic heat conduction equation into four real first-order differential equations, followed by employment of the Bulirsch-Stoer method with adoptive stepsize. The applicability of the procedure is then demonstrated in few issues, important for the understanding of thermal signal distribution.