<p>In complex physical systems, conventional differential equations fall short in capturing non-local and memory effects. Fractional differential equations (FDEs) effectively model long-range interactions with fewer parameters. However, deriving FDEs from physical principles remains a significant challenge. This study introduces a stepwise data-driven framework to discover explicit expressions of FDEs directly from data. The proposed framework combines deep neural networks for data reconstruction and automatic differentiation with Gauss-Jacobi quadrature for fractional derivative approximation, effectively handling singularities while achieving fast, high-precision computations across large temporal/spatial scales. To optimize both linear coefficients and the nonlinear fractional orders, we employ an alternating optimization approach that combines sparse regression with global optimization techniques. We validate the framework on various datasets, including synthetic anomalous diffusion data, experimental data on the creep behavior of frozen soils, and single-particle trajectories modeled by Lévy motion. Results demonstrate the framework’s robustness in identifying FDE structures across diverse noise levels and its ability to capture integer-order dynamics, offering a flexible approach for modeling memory effects in complex systems.</p>