In computational mathematics, when dealing with a large linear discrete problem (e.g., a linear system) arising from the numericaldiscretization of a differential equation (DE), knowledge of the spectral distribution of the associated matrix has proved to be useful information for designing/analyzing appropriate solvers–-especially, preconditioned Krylov and multigrid solvers–-for the considered problem. Actually, this spectral information is of interest also in itself as long as the eigenvalues of the aforementioned matrix represent physical quantities of interest, which is the case for several problems from engineering and applied sciences (e.g., the study of natural vibration frequencies in an elastic material).The theory of generalized locally Toeplitz (GLT) sequences is a powerful apparatus for computing the asymptotic spectral distribution of matrices $A_n$ arising from virtually any kind of numerical discretization of DEs. Indeed, when the mesh-fineness parameter $n$ tends to infinity, these matrices $A_n$ give rise to a sequence $\\{A_n\\}_n$, which often turns out to be a GLT sequence or one of its “relatives”, i.e., a block GLT sequence or a reduced GLT sequence. In particular, block GLT sequences are encountered in the discretization of systems of DEs as well as in the higher-order finite element or discontinuous Galerkin approximation of scalar/vectorial DEs.This work is a review, refinement, extension, and systematic exposition of the theory of block GLT sequences. It also includes several emblematic applications of this theory in the context of DE discretizations.