For diffraction-limited optical systems, an accurate physical optics model is necessary to properly evaluate instrument performance. Astronomical observatories outfitted with coronagraphs for direct exoplanet imaging require physical optics models to simulate the effects of misalignment and diffraction. Accurate knowledge of the observatory’s point-spread function (PSF) is integral for the design of high-contrast imaging instruments and simulation of astrophysical observations. The state of the art is to model the misalignment, ray aberration, and diffraction across multiple software packages, which complicates the design process. Gaussian beamlet decomposition (GBD) is a ray-based method of diffraction calculation that has been widely implemented in commercial optical design software. By performing the coherent calculation with data from the ray model of the observatory, the ray aberration errors can be fed directly into the physical optics model of the coronagraph, enabling a more integrated model of the observatory. We develop a formal algorithm for the transfer-matrix method of GBD and evaluate it against analytical results and a traditional physical optics model to assess the suitability of GBD for high-contrast imaging simulations. Our GBD simulations of the observatory PSF, when compared to the analytical Airy function, have a sum-normalized RMS difference of ≈10 ^{− 6}. These fields are then propagated through a Fraunhofer model of an exoplanet imaging coronagraph where the mean residual numerical contrast is 4 × ^{10 − 11}, with a maximum near the inner working angle at 5 × 10 ^{− 9}. These results show considerable promise for the future development of GBD as a viable propagation technique in high-contrast imaging. We developed this algorithm in an open-source software package and outlined a path for its continued development to increase the accuracy and flexibility of diffraction simulations using GBD.