This paper investigates the numerical solution of a three-dimensional time–space fractional functional partial differential equation involving Caputo fractional derivatives in time and Riesz fractional derivatives in space. Owing to the inherent difficulties associated with fractional operators, such as weakly singular kernels and strong nonlocality, conventional numerical methods often encounter limitations in terms of accuracy and computational efficiency. To address these challenges, a hybrid numerical scheme is developed in which the Caputo fractional derivative is discretized by means of cubic spline interpolation, providing high-order tem-poral accuracy, while advanced meshless techniques are employed for the spatial approximation of the Riesz fractional derivative, thereby effectively managing its global support and singular behavior. The resulting fully discrete scheme is rigorously analyzed, and its unconditional stability is established using the energy method. Extensive numerical experiments on various three-dimensional domains demonstrate the high ac-curacy of the proposed method, with convergence rates consistent with the theoretical analysis, as well as favorable computational efficiency. Overall, the proposed approach offers a robust and flexible computa-tional framework for the numerical treatment of complex fractional models arising in physics, engineering, and applied sciences.
M.J. Huntul (Mon,) studied this question.