We develop a novel and efficient boundary integral equation method based on the spatio-temporal formulation for the two-dimensional dynamic and quasistatic analyses of an earthquake fault in a single scheme. A major advantage of this method is its applicability to the analysis of non-planar faults with the same degree of accuracy as to that of planar faults. Calculation time and memory requirement are reduced through the employment of asymptotic representations of the integration kernels appearing in the convolution integral. Asymptotic kernels are factorized into terms dependent on space or time alone, resulting in efficient numerical computations. In addition, the dependence on time is found to vanish in the asymptotic kernels far behind the S-wave front, which also contributes to the time-saving efficiency of the calculations. We show that, in a dynamic analysis, if a 3% error is allowed for the slip rate, computation time and memory requirement are reduced by 25% and 45%, respectively, in an in-plane fault case, and by 60% and 75%, respectively, in an anti-plane fault case. This method can be employed as a powerful numerical tool in simulating an entire earthquake cycle consisting of both quasi-static and dynamic processes with a more realistic non-planar geometry of faults.