The design and optimization of medical ultrasound modalities require a method for modeling the propagation of pressure wave fields through biomedical tissue. This method needs to be capable of modeling various phenomena, e.g. nonlinear propagation, attenuation, and scattering caused by arbitrary inhomogeneities in the acoustic properties of the medium. The Iterative Nonlinear Contrast Source (INCS) method is a full wave method which makes use of a contrast source approach to model nonlinear propagation of ultrasound. Originally, the method was used to model nonlinear propagation in homogeneous nonlinear media with frequency power law attenuation. In this thesis, the basic steps of INCS are conserved and the method is extended to deal with spatially varying attenuation, coefficient of nonlinearity and speed of sound. The method developed uses a generalized form of the Westervelt equation which includes a spatially dependent coefficient of nonlinearity, speed of sound and relaxation function. The generalized Westervelt equation is recast into an integral equation and solved using three different iterative schemes: Neumann, Bi-CGSTAB and Steepest Descent. The performances of these schemes are analyzed and compared. Furthermore, as an example of a possible application of the method developed, a feasibility study of a new beam forming technique, i.e. parallel transmit beam forming using orthogonal frequency division multiplexing applied to harmonic imaging, is performed. Numerical studies show the capability of this technique to reduce the presence of unwanted side-lobes and meanwhile increase the amplitude of the main beam as compared to standard parallel beam forming in reception. These improvements are expected to positively influence the signal to noise ratio and the achievable penetration depth of a given imaging system. A drawback of the technique proposed is a reduction of the axial resolution due to utilization of pulses with a narrower bandwidth. Finally, measurements in water are performed which confirmed the feasibility of the technique proposed for a practical transducer.