This paper develops a high-order, fully meshless numerical framework for solving strongly coupled nonlinear two-dimensional reaction–diffusion systems governed by distributed-order fractional derivatives in both time (Caputo type) and space (Riesz type), with additional nonlocal integral interactions. A linearly implicit L2- \(1_\sigma \) –Gauss–Legendre time semi-discretization is first introduced that achieves a uniform convergence rate of \(3-\alpha _{\max } \in [2,3)\) , where \(\alpha _{\max }:=\max _{1\le k\le K}\alpha _k,\) and \(\alpha _k\) denotes the fractional-order quadrature node associated with the distributed-order Caputo fractional derivative. The distributed-order Riesz spatial operators are subsequently discretised in a completely meshless fashion by combining Gauss–Legendre quadrature in the spatial order variable \(\beta \in [0,2]\) with high-order polyharmonic-spline RBF-FD stencils, yielding spectral accuracy in \(\beta \) and algebraic order \(d+1\) on arbitrary scattered nodes. The resulting fully discrete scheme is proved to be unconditionally stable and optimally convergent with the sharp rate \(O(\tau ^{3-\alpha _{\max }} + h^{d+1} + 4^{-M} + 4^{-K}),\) where \(\tau \) denotes the temporal step size, h is the spatial fill distance associated with the scattered nodes, and M and K represent the Gauss–Legendre quadrature orders in the distributed spatial-order and temporal-order variables, respectively, in the discrete \(\ell _h^2\) -norm using a novel energy argument that simultaneously controls the distributed-order memory, nonlocal terms, and globally Lipschitz nonlinear reactions. The theoretical findings are rigorously validated by two representative examples: (i) a manufactured smooth solution with known closed-form expression that confirms the predicted convergence orders, and (ii) a physically motivated nonlinear system with highly singular weight \(\Gamma (7/2-\alpha )\) and nonlocal coupling, for which no analytical solution exists, demonstrating the scheme exhibits excellent robustness and accurately resolves complex spatiotemporal patterns dominated by strong memory effects. The proposed method therefore provides a powerful, flexible, and uniformly accurate tool for the numerical investigation of complex distributed-order fractional PDEs arising in anomalous transport, viscoelasticity, and nonlocal interaction phenomena.
