In this article, we investigate the utilization of Riemann–Liouville’s fractional integral and the Caputo derivative in the application of the Optimal Auxiliary Function Method (OAFM). The extended OAFM is employed to analyze fractional non-linear coupled ITO systems and non-linear KDV systems, which feature equations of a fractional order in time. We compare the results obtained for the ITO system with those derived from the Homotopy Perturbation Method (HPM) and the New Iterative Method (NIM), and for the KDV system with the Laplace Adomian Decomposition Method (LADM). OAFM demonstrates remarkable convergence with a single iteration, rendering it highly effective. In contrast to other existing analytical approaches, OAFM emerges as a dependable and efficient methodology, delivering high-precision solutions for intricate problems while saving both computational resources and time. Our results indicate superior accuracy with OAFM in comparison to HPM, NIM, and LADM. Additionally, we enhance the accuracy of OAFM through the introduction of supplementary auxiliary functions.