In this study, we develop a continuous Galerkin approximation method for solving a class of linear third-kind Volterra integral equations (VIEs), which arise in modeling physical and biological phenomena and exhibit singular behavior due to the vanishing coefficient \(t^{\beta }\) . Our approach employs the shifted Legendre polynomials as basis functions within a piecewise polynomial space, combining their orthogonality properties and spectral accuracy to handle the inherent singularities efficiently. The method transforms the integral equation into an algebraic system via variational formulation, ensuring solution continuity across subintervals while maintaining computational stability. We provide a rigorous convergence analysis, proving that the scheme achieves optimal convergence rates of \(\mathcal {O}\left( h^{m+1}\right) \) for odd-degree polynomials and \(\mathcal {O}\left( h^{m}\right) \) for even degrees, where m is the polynomial order and h the mesh size, provided the exact solution is sufficiently smooth. For non-smooth solutions, the method automatically adapts to achieve the best possible algebraic convergence rate permitted by the solution’s regularity. Numerical experiments validate the theoretical results, demonstrating high-order accuracy for both smooth trigonometric and exponential solutions, while also confirming the method’s robustness for non-smooth cases where convergence becomes regularity-limited. This analysis extends prior work on third-kind VIEs by comprehensively addressing the challenges posed by both smooth and non-smooth kernels and vanishing coefficients.