The motivation of this work is to analyze the in-band intermodulation distortion (IMD) occurring in surface acoustic wave (SAW) devices, using a recently developed fast method based on the input-output equivalent sources (IOES). The method calculates the equivalent current sources of a given harmonic (H) or IMD, which when applied at the boundaries of any uniform nonlinear region produce the same nonlinearities as the full distributed circuit. The accuracy of the method is validated with a very simplified SAW resonator with ten digits, which is modeled by a discretized Mason-based circuit. The IOES method provides equal results to the ones obtained through harmonic balance (HB) simulations, performed by means of commercial software, being the first 1000 times faster. Once the accuracy of the method is guaranteed, it is used to analyze the measured in-band IMD3 of several lithium tantalite 42° cut leaky SAW (LSAW) resonators with different pitches and duty factors at the B66 long term evolution (LTE) frequency band. Those resonators are comprised of 100 and 20 electrode pairs for the active region and each of the reflectors, respectively, which implies the analysis of a very large distributed nonlinear problem with thousands of nonlinear local sources. The IOES method takes 35.4 s in simulating 51 frequency points, whereas this simulation is not possible using a commercial HB simulator on a general-purpose computer.