In nondestructive testing, the generation of higher harmonics and the mixing of elastic waves are used to measure the nonlinearity parameters, which in turn are closely related to the microscopic state of the material. A 3-D numerical tool that can simulate large scale, four-dimensional, nonlinear elastic wavefields would be very useful for the reliable interpretation of experimental results. Here, a method is presented that can efficiently compute the nonlinear displacement fields in a homogeneous, isotropic elastic medium, taking into account its third-order elastic constants (A, B, C). The method is based on the Neumann iterative solution of an integral equation involving a Green's function of the linear `background' medium, and a contrast source representing the nonlinearity of the medium. The integral equation is solved iteratively, and the computations are based on Fast Fourier Transforms using a sampling rate close to the Nyquist limit, i.e., two grid points of the shortest wavelength. The displacement fields are evaluated using the scalar and vector potential functions representing the compressional and shear displacements. In this presentation, the suitability of the INCS method for modeling the harmonics of 3-D nonlinear elastic waves and its validation by comparison with an analytical benchmark solution, will be presented.