There are many applications in which the properties of nonlinear waves play a role. For example, in medical diagnostic imaging, the higher harmonics in nonlinear ultrasound are used to reduce clutter from nearby artifacts and to improve resolution, and in non-destructive testing nonlinear mixing products from crossing ultrasound beams are used to assess the state of materials. Simulation of nonlinear waves is important for optimizing equipment and quantitatively assessing phenomena. Unfortunately, analytical approaches cannot deal with complex realistic situations, while traditional Finite Element and Finite Difference methods require extremely large grids to capture waves with many higher harmonics in large three-dimensional domains. The Iterative Nonlinear Contrast Source method was developed to overcome these issues. It is based on an integral equation involving the analytic Green’s function of the linearized medium and a distributed contrast source representing nonlinearity. The integral equation is solved iteratively, and the computations are based on Fast Fourier Transforms that only require two grid points of the shortest wavelength. In this presentation, the basic principles behind the method will be explained and its extension to lossy and inhomogeneous media will be shown. Moreover, a novel extension will be presented that enables the computation of nonlinear elastic waves.