This paper presents an efficient and accurate numerical approach for solving fractional differential equations (FDEs), particularly those with delay terms, using a truncated series of Fermat polynomials. The developed method converts the fractional differential equation and its associated initial conditions into a system of algebraic equations by employing the Galerkin spectral technique combined with the operational matrix of fractional-order derivatives in the Caputo sense. Fermat polynomials are adopted as basis functions due to their advantageous analytical properties and recursive structure, making them suitable for constructing spectral approximations. The operational matrix formulation significantly reduces computational complexity and facilitates efficient implementation. The approximate solution is obtained with high accuracy using only a limited number of basis functions, ensuring rapid convergence. To demonstrate the effectiveness of this approach, several numerical examples, including fractional delay differential equations, are investigated. The results show excellent agreement with exact solutions, even for non-polynomial problems, and comparisons with existing methods confirm the superiority, stability, and computational efficiency of the proposed technique.