Manipulating the incident wavefront in biomedical applications to enhance the penetration depth and energy delivery in scattering media such as biological tissue has gained a lot of attention in recent years. However, focusing inside scattering media and examining the electromagnetic field inside the medium still is an elaborate task. This is where electromagnetic field simulations that model the wavefront shaping process can help us understand how the focal near field evolves at different depths. Here we use a two-step beam synthesis method to simulate the scattering of complex incident wavefronts by well-characterized media. The approach uses plane wave electromagnetic near-field solutions in combination with an angular spectrum approach to model different light beams. We apply this approach to various two-dimensional scattering media and investigate the focus intensity over depth while scanning with and without phase optimization. We find that the scanned non-optimized beams have two regions characterized by exponential decays. The absolute progression of the focus intensity over depth for phase-optimized beams using all channels can be described by solutions of the radiative transfer theory. Furthermore, the average enhancement factor over depth of the phase-optimized focus intensity compared to that without optimization is investigated for different numerical apertures and scattering media. Our results show that, albeit the incident beam is diffusively scattered, the theoretical enhancement for a large number of optimization channels cannot be reached due to correlations between the channels. An increase in focus depth and an increase in the numerical aperture reduces the difference between the expected theoretical and simulated enhancement factors.