We propose a novel method for the wave-optical simulation of diffractive optical elements (DOEs) like metasurfaces or computer-generated holograms (CGHs). Existing techniques mostly rely on the assumption of local periodicity to predict the performance of elements. The utilization of a specially adapted finite-difference beam propagation method (BPM) allows the semi-rigorous simulation of entire DOEs within a reasonable runtime due to linear scaling with the number of grid points. Its applicability is demonstrated by the simulation of a metalens and a polarization-dependent beamsplitter, both based on effective-medium metasurfaces. A comparison shows high conformity to rigorous simulations.