Improvements in industrial seismic, seismological, acoustic, or interferometric theory and applications often result in quite subtle changes in sound quality, seismic images, or information which are nevertheless crucial for improved interpretation or experience. When evaluating new theories and algorithms using synthetic data, an important aspect of related research is therefore that numerical errors due to wavefield modeling are reduced to a minimum. We present a new MATLAB code based on the Foldy method that models theoretically exact, direct, and scattered parts of a wavefield. Its main advantage lies in the fact that while all multiple scattering interactions are taken into account, unlike finite-difference or finite-element methods, numerical dispersion errors are avoided. The method is therefore ideal for testing new theory in industrial seismics, seismology, acoustics, and in wavefield interferometry in particular because the latter is particularly sensitive to the dynamics of scattering interactions. We present the theory behind the Foldy acoustic modeling method and provide examples of its implementation. We also benchmark the code against a good finite-difference code. Because our Foldy code was written and optimized to test new theory in seismic interferometry, examples of its application to seismic interferometry are also presented, showing its validity and importance when exact modeling results are needed.