Asynchronous vibration was generated between the main bridge and approach spans or abutments due to differences in stiffness and mass during an earthquake, thus further leading to pounding at the bilateral beam ends. By taking a T-shaped rigid frame bridge as an example, the bilateral pounding model was abstracted, and the earthquake response spectra considering pounding at the bilateral beam ends were studied, including the maximum displacement spectrum, the acceleration dynamic coefficient spectrum, the pounding force response spectrum, and the response spectrum for the number of pounding events. An improved precise pounding algorithm was proposed to solve the dynamic equation of the bilateral pounding model. This algorithm is based on the precise integration method for solving the second-order dynamic differential equation and reduces the order thereof by introducing a new velocity vector and uses the series method to find the nonhomogeneous term. The system matrix is simpler, and the inversion of the system matrix can be avoided. On this basis, a multipoint earthquake-induced pounding response spectrum program was developed. A total of 18 seismic waves from Class II sites were selected, and the response spectra of 18 waves were analyzed using this new program. Furthermore, the effects of structural stiffness, mass, stiffness of contact element, pounding recovery coefficient, and peak ground acceleration (PGA) on the earthquake response spectrum were studied. Through the analysis of earthquake response spectra and a parametric study, the phenomenon of earthquake-induced pounding of bridges was clarified to the benefit of the analysis and engineering control of earthquake-induced pounding of bridges.