Not sure what your application is, but you may not need POP if diffraction effects can be neglected (apart from propagation from the exit pupil to the image, or focal, plane -- which can be handled with the FFT and/or Huygens PSF analysis).
However, if you do need to account for high-NA diffraction (e.g., you implement any form of spatial filtering), or you simply need to know the wavefront profile at various locations along the propagation path, then you may want to take a look at the following software package: DIFFRACT, authored by Prof. Masud Mansuripur, who is one of the world’s foremost experts on physical optics. The software was originally written to support simulation of high-NA optical disk drives. The interface is a bit dated now, but the numerical computations are very accurate. It provides three user-selectable propagation regimes (near-field, far-field, and the intermediate range), none of which rely on the paraxial Fresnel approximation. Thin optical elements can be inserted in the path (e.g., thin lenses with optional user-supplied Zernike or Seidel aberration, phase plates, apertures, gratings, DOEs, etc.), and ray tracing with thick elements is also supported.
The mathematical details are published in the following papers:
Certain computational aspects of vector diffraction problems
Distribution of light at and near the focus of high-numerical-aperture objectives
Errata