This chapter provides a comprehensive presentation of global numerical methods for solving FODEs employing the polynomial and non-polynomial bases, introduced in Chapter 2. The FODEs of interest will be initial-/boundary-value problems, posed using a variety of fractional derivatives (e.g., Caputo, Riemann–Liouville, Riesz, one-sided, two-sided, variable-order, distributed order, etc.), introduced in Chapters 1 and 2. We devote Sections 3.1 and 3.2 to introducing a series of variational and non-variational spectral methods in single domains, where the solution singularities can occur at the initial or boundary points. In a variational formulation of an FODE, one first obtains the weak (variational) form of the given equation, where the highest derivative order is reduced using integration-by-parts, and then solves the variational formulation by constructing the corresponding (finite-dimensional) solution and test subspaces. In non-variational problems, one rather directly solves the strong (original) FODE, hence assuming a higher regularity in the solution. Moreover, we introduce spectral element methods (SEM) for FODEs in multiple domains for the main purpose of capturing possible interior/boundary singularities.