In this paper we present a robust pseudo-random method for propagating an electromagnetic wavefront through an arbitrary optical system. The wavefront at an arbitrary plane is obtained by the discrete sampling of the wavefront in regular regions equally distributed on the entrance pupil and the appropriate modelling of the optical properties of the system. The discretization permits us to treat each region as a plane wave, so long as the area is small compared to the area of the pupil, therefore allowing us to apply the electromagnetic approximations of refraction and reflection during the transfer through an optical system. We can therefore account for amplitude and phase modulation of the wavefront due to the optical system, without making any assumptions about the shape of the optical elements. Furthermore, our numerical integration method on an arbitrary plane avoids singularities due to the classical analytical integrals, while still obtaining results comparable to rigorous electromagnetic theory. We have applied the method to simulating the propagation of both plane waves and spherical waves. The well known interference patterns of classical experiments such as Young's interference fringes or Newton's rings were reproduced accurately, with respect to results obtained applying analytical methods. We then successfully applied the method to analyze a Michelson interferometer set-up, demonstrating the robustness of the calculations. Since the propagation of the wavefront is possible with this method, in the future we plan to apply the method to simulating electrically large diffractive optical elements within a complex optical system, for which rigorous analytical methods may not be available, and other numerical methods generally require large computer resources.