A time domain boundary integral equation method, which is named as quasi-initial condition method, is applied in this paper to solve the transient heat conduction problem. In conventional implementations, however, this method suffers from a numerically unstable problem when the time step is small. To improve the numerical stability of the method, a time step amplification method is proposed. In the proposed method, an amplified time step is adopted to compute the temperature and the flux at the virtual time point. The boundary condition at that virtual time point is determined through a linear interpolation by the conditions at the current time step point and the quasi-initial time. Furthermore, the heat generation in the virtual time step is assumed to be constant which is the same as that in the real time step. The temperature and the flux at the current step time point are then computed through a linear interpolation over the time interval. A short but not rigorous deduction of this method is presented to show that this method is valid in solution to problems in which the temperature and the flux vary linearly respect to time. Numerical examples further demonstrate the numerical stability of the proposed method. (C) 2015 Elsevier Ltd. All rights reserved.