A novel, free from paraxial approximation and computationally efficient numerical algorithm capable of predicting 4D acoustic fields in lossy and nonlinear media from arbitrary shaped sources (relevant to probes used in medical ultrasonic imaging and therapeutic systems) is described. The new WE (wave envelopes) approach to nonlinear propagation modeling is based on the solution of the second order nonlinear differential wave equation reported in [J. Wójcik, J. Acoust. Soc. Am. 104 (1998) 2654–2663; V.P. Kuznetsov, Akust. Zh. 16 (1970) 548–553]. An incremental stepping scheme allows for forward wave propagation. The operator-splitting method accounts independently for the effects of full diffraction, absorption and nonlinear interactions of harmonics. The WE method represents the propagating pulsed acoustic wave as a superposition of wavelet-like sinusoidal pulses with carrier frequencies being the harmonics of the boundary tone burst disturbance. The model is valid for lossy media, arbitrarily shaped plane and focused sources, accounts for the effects of diffraction and can be applied to continuous as well as to pulsed waves. Depending on the source geometry, level of nonlinearity and frequency bandwidth, in comparison with the conventional approach the Time-Averaged Wave Envelopes (TAWE) method shortens computational time of the full 4D nonlinear field calculation by at least an order of magnitude; thus, predictions of nonlinear beam propagation from complex sources (such as phased arrays) can be available within 30–60 min using only a standard PC. The approximate ratio between the computational time costs obtained by using the TAWE method and the conventional approach in calculations of the nonlinear interactions is proportional to 1/N∼2, and in memory consumption to 1/N∼ where N∼ is the average bandwidth of the individual wavelets. Numerical computations comparing the spatial field distributions obtained by using both the TAWE method and the conventional approach (based on a Fourier series representation of the propagating wave) are given for circular source geometry, which represents the most challenging case from the computational time point of view. For two cases, short (2 cycle) and long (8 cycle) 2 MHz bursts, the computational times were 10 min and 15 min versus 2 h and 8 h for the TAWE method versus the conventional method, respectively.