When simulating phenomena in physics, engineering, or applied sciences, often one has to deal with functional equations that do not admit an analytical solution. Describing these real situations is, however, possible, resorting to one of its numerical approximations and treating the resulting mathematical representation. This thesis is placed in this context: Indeed the purpose is that of furnishing several useful tools to deal with some computational problems, stemming from discretization techniques. In most of the cases the numerical methods we analyse are the classical $\mathbb{Q}_p$ Lagrangian FEM and the more recent Galerkin B-spline Isogeometric Analysis (IgA) approximation and Staggered Discontinuous Galerkin (DG) methods. As our model PDE, we consider classical second-order elliptic differential equations and the Incompressible Navier-Stokes equations. In all these situations the resulting matrix sequences $\{A_n\}_{n}$ possess a structure, namely they belong to the class of Toeplitz matrix sequences or to the more general class of Generalized Locally Toeplitz (GLT) matrix sequences, in the most general block $k$-level case. Consequently, the spectral analysis of the coefficient matrices plays a crucial role for an efficient and fast resolution. Indeed the convergence properties of iterative methods proposed, like multigrid or preconditioned Krylov techniques, are strictly related to the notion of symbol of the coefficient matrix sequence. In our setting the symbol is a function which asymptotically provides a reasonable approximation of the eigenvalues [singular values] of $A_n$ by its evaluations of an uniform grid on its domain. These reasons, and many others, make the research of more and more efficient eigensolvers relevant and topical. In this direction, the second goal of this thesis is to provide new tools for computing the spectrum of preconditioned banded symmetric Toeplitz matrices, Toeplitz-like matrices, $n^{-1}K_n^{[p]}$, $nM_n^{[p]}$, $n^{-2}L_n^{[p]}$, coming from the B-spline IgA approximation of $-u"= \lambda u$, plus its multivariate counterpart for {$-\Delta{u}= \lambda u$}, and block and preconditioned block banded symmetric Toeplitz matrices. For all the above cases we propose new algorithms based on the classical concept of symbol, but with an innovative view on the errors of the approximation of eigenvalues by the uniform sampling of the symbol. The algorithms devised are special interpolation-extrapolation procedures performed with a high level of accuracy and only at the cost of computing of the eigenvalues of a moderate number of small sized matrices.