End-to-end simulation of the influence of the optical train on the observed scene is important across optics and is particularly important for predicting the science yield of astronomical telescopes. As a consequence of their goal of suppressing starlight, coronagraphic instruments for high-contrast imaging have particularly complex field-dependent point-spread-functions (PSFs). The Roman Coronagraph Instrument (CGI), Hybrid Lyot Coronagraph (HLC) is one example. The purpose of the HLC is to image exoplanets and exozodiacal dust in order to understand dynamics of solar systems. This paper details how images of exoplanets and exozodiacal dust are simulated using some of the most recent PSFs generated for the CGI HLC imaging mode. First, PSFs are generated using physical optics propagation techniques. Then, the angular offset of pixels in image scenes, such as exozodiacal dust models, are used to create a library of interpolated PSFs using interpolation and rotation techniques, such that the interpolated PSFs correspond to angular offsets of the pixels. This means interpolation needs only be done once and an image can then be simulated by multiplying the vector array of the model astrophysical scene by the matrix array of the interpolated PSF data. This substantially reduces the time required to generate image simulations by reducing the process to matrix multiplication, allowing for faster scene analysis. We will detail the steps required to generate coronagraphic scenes, quantify the speed-up of our matrix approach versus other implementations, and provide example code for users who wish to simulate their own scenes using publicly available HLC PSFs.