The problem of heat conduction in composite media is formulated in terms of integral equations by use of the tempered elementary solution of the unsteady diffusion operator. The integral equations are discretized in both time and space. Each homogeneous subre-gion and its surface and interfaces are represented by tetrahedral and triangular finite elements with quadratic variation of geometry and of unknown functions with respect to intrinsic coordinates. The time dependence of the unknown functions is assumed to be linear over each time step. Integration in the time domain is performed analytically, while integration with respect to spatial coordinates is developed by numerical quadrature formulas. As a result, a linear algebraic system of equations for the unknown in the nodal points is obtained. The procedure was checked with a number of test cases. The numerical experiments show that good accuracy may be obtained provided the parameters of the computation, in particular the number of Gaussian quadrature points and the time step for a given nodal mesh, are suitably chosen. Proper choice of these parameters is particularly important when dealing with composite media in which the thermal diffusivtiy varies widely.