Numerical methods capable of modeling nonlinear pressure wave fields propagating through inhomogeneous biomedical tissue are essential for the design and optimization of ultrasound transducers or devices. The Iterative Nonlinear Contrast Source (INCS) method is an accurate method for modeling three-dimensional nonlinear acoustic wave fields. Originally it was capable of modeling nonlinear wave fields in homogeneous lossy tissue. Recently, it has been extended to deal with spatially varying coefficient of nonlinearity and attenuation. The method recasts a generalized form of the Westervelt equation into an integral equation which was originally solved using a Neumann scheme. This scheme allows to model moderate losses and non-linearity. Problems with the convergence may occur for realistic speed of sound contrast. Here, we present a different solution method which makes it possible to treat, besides spatially varying coefficient of nonlinearity and attenuation, also realistic speed of sound contrasts. The nonlinear integral equation is solved using a steepest descent scheme. The method has been used to compute the three-dimensional nonlinear pressure wave field generated by a 40 element linear array and propagating through a medium with spatially varying coefficient of nonlinearity, attenuation and speed of sound. Simulations have been performed up to the 5th harmonic component.