Fiber-optic probes for sensing and biomedical imaging applications such as optical coherence tomography (OCT) frequently employ sections of graded-index (GRIN) fiber to re-focus the diverging light from the delivery fiber. Such GRIN fiber microlenses often possess aberrations that cause significant distortions of the focused output beam. Current design methods based on ABCD matrix transformations of Gaussian beams cannot model such effects and are therefore inadequate for the analysis and design of high-performance probes that require diffraction-limited output beams. We demonstrate use of the beam propagation method (BPM) to analyze beam distortion in GRIN-lensed fibers resulting from index profiles that exhibit a deviation from the ideal parabolic shape or artifacts such as ripples or a central dip. Furthermore, we demonstrate the power of this method for exploring novel probe designs that incorporate GRIN phase masks to generate wavefront-shaped output beams with extended depth-of-focus (DOF). We present results using our method that are in good agreement with experimental data. The BPM enables accurate simulation of fiber probes using non-ideal or custom-engineered GRIN fibers with arbitrary refractive index profiles, which is important in the design of high-performance fiber-based micro-imaging systems for biomedical applications.