Many novel medical echography modalities, like super harmonic imaging (SHI), employ imaging of higher harmonics arising from nonlinear propagation. Design optimization of dedicated imaging probes for these modalities requires accurate computation of the higher harmonic beam profiles. The existing iterative nonlinear contrast source method can perform such simulations. With this method, the nonlinear term of the Westervelt equation is considered to represent a distributed contrast source in a linear background medium. The full transient nonlinear acoustic wave field follows from the Neumann iterative solution of the corresponding integral equation. Each iteration involves the spatiotemporal convolution of the background Green's function with an estimate of the contrast source, and appropriate filtering enables a discretization approaching two points per smallest wavelength or period. To further reduce the computational load and to anticipate extension of the method, in this presentation, the effect of linearizing the nonlinear contrast source is investigated. This enables the replacement of the Neumann scheme by more efficient schemes. Numerical simulations show that for design purposes the effect of linearization on the harmonic components involved with SHI remains sufficiently small. Moreover, it is shown that a significant decrease in computational costs may be achieved by using a Bi-CGSTAB scheme.