This work presents a new method to analyze weak distributed nonlinear (NL) effects, with a focus on the generation of harmonics (H) and intermodulation products (IMD) in bulk acoustic wave (BAW) resonators and filters composed of them. The method consists of finding equivalent current sources [input-output equivalent sources (IOES)] at the H or IMD frequencies of interest that are applied to the boundary nodes of any layer that can contribute to the nonlinearities according to its local NL constitutive equations. The new methodology is compared with the harmonic balance (HB) analysis, by means of a commercial tool, of a discretized NL Mason model, which is the most used model for NL BAW resonators. While the computation time is drastically reduced, the results are fully identical. For the simulation of a seventh-order filter, the IOES method is around 700 times faster than the HB simulations.