This paper proposes a novel and efficient spectral numerical framework for approximating solutions of Volterra–Fredholm fractional integro-differential equations with variable-order derivatives. The method is based on the use of normalized and shifted Vieta–Lucas orthogonal polynomials as basis functions to represent the unknown solution. New operational matrices are derived for both the Caputo fractional derivative and the Riemann–Liouville fractional integral operators, enabling the transformation of the original fractional system into a system of linear or nonlinear algebraic equations. The proposed scheme exhibits favorable properties such as simple implementation, strong numerical stability, and reliable convergence behavior. Rigorous theoretical results concerning the existence, uniqueness, and convergence of the approximate solution are established. Several numerical examples are presented to demonstrate the accuracy and computational efficiency of the proposed approach.
Mizani et al. (Sun,) studied this question.