In this work, we develop a Jacobi spectral Galerkin method for solving weakly singular systems of second-kind Volterra integro-differential equations. To achieve high-order accuracy, suitable variable and function transformations are employed to map the original problem onto the standard interval \([-1, t]\) . Rigorous existence and uniqueness results for the solutions of systems of second-kind Volterra integro-differential equations are also established. Alongside the theoretical framework, we present a detailed numerical algorithm for practical implementation. The convergence of the proposed method is analyzed thoroughly, with error estimates derived both in the infinity norm and in a weighted \(L^2\) -norm. A distinct feature of this work is the comprehensive convergence analysis, which sets it apart from existing studies and highlights the robustness and precision of the method. Finally, the theoretical predictions are validated through numerical experiments.
