An efficient numerical solution is proposed to solve transient heat conduction problems in three dimensions using the boundary element method without volume discretization. The thermal conductivity is assumed to be constant, and no internal heat sources or sinks are present. The boundary integral equation with time-dependent fundamental solution of the governing equation is discretized by using constant time steps. Then, the time integrations of the integral kernels result in incomplete gamma functions (IGF). Using truncated series expansion of the IGFs, the time and space variables are separated in multiplicative factors. Thus, the product of all the spatial quantities in the BIE can be integrated once and saved during the computation for each term of the series expansion. The matrices at each time step can be generated from the stored spatial quantities by multiplication by time factors. The solution scheme presented allows for dynamically changing the step as the solution evolves. The question of disk storage and spatial integrations are discussed. Three numerical examples are presented to demonstrate the efficiency of the present scheme and substantial reduction of the CPU requirements.