A numerical method is developed to deconvolve complex body waves into a multiple shock sequence. The method assumes that all constituent events have identical fault geometry and depth, and the far-field source time function is obtained as a superposition of ramp functions. The height and onset time of the ramp functions are determined by matching synthetic and observed waveforms in the least-square sense. The individual events are identified by pairs of ramp functions or discrete trapezoidal pulses in the source time sequence. The method is applicable to both single and multi-station data. The 1976 Guatemala earthquake is analyzed as a test case, demonstrating the method's effectiveness in systematically analyzing multiple event sequences. The method involves modeling the far-field source time function as a convolution of dislocation velocity and fault area expansion rate. The source time function is expressed as a superposition of ramp functions, with each event represented by a pulse. The method is applied to single-station and multi-station data, with the latter requiring additional parameters for source location. The method is validated using the Guatemala earthquake data, where the source time function is decomposed into multiple pulses, and the seismic moment of individual events is estimated. The results show that the method can effectively identify multiple events and their characteristics, including the timing and location of ruptures. The method is useful for analyzing complex earthquake sequences and provides insights into fault behavior and rupture dynamics.A numerical method is developed to deconvolve complex body waves into a multiple shock sequence. The method assumes that all constituent events have identical fault geometry and depth, and the far-field source time function is obtained as a superposition of ramp functions. The height and onset time of the ramp functions are determined by matching synthetic and observed waveforms in the least-square sense. The individual events are identified by pairs of ramp functions or discrete trapezoidal pulses in the source time sequence. The method is applicable to both single and multi-station data. The 1976 Guatemala earthquake is analyzed as a test case, demonstrating the method's effectiveness in systematically analyzing multiple event sequences. The method involves modeling the far-field source time function as a convolution of dislocation velocity and fault area expansion rate. The source time function is expressed as a superposition of ramp functions, with each event represented by a pulse. The method is applied to single-station and multi-station data, with the latter requiring additional parameters for source location. The method is validated using the Guatemala earthquake data, where the source time function is decomposed into multiple pulses, and the seismic moment of individual events is estimated. The results show that the method can effectively identify multiple events and their characteristics, including the timing and location of ruptures. The method is useful for analyzing complex earthquake sequences and provides insights into fault behavior and rupture dynamics.