In this paper, a new approach for numerically solving the system of fractional integrodifferential equations is devised. To approximate the issue, we employ Vieta–Fibonacci polynomials as basis functions and derive the projection method for Caputo fractional order for the first time. An efficient transformation reduces the problem to a system of two independent equations. Solving two algebraic equations yields an approximate solution to the problem. The proposed method’s efficiency and accuracy are validated. We demonstrate the existence of the solution to the approximate problem and conduct an error analysis. Numerical tests reinforce the interpretations of the theory.