We introduce a boundary penalization technique to improve the spectral approximation of isogeometric analysis (IGA). The technique removes the outliers appearing in the highfrequency region of the approximate spectrum when using the C p−1 , p-th (p ≥ 3) order isogeometric elements. We focus on the classical Laplacian (Dirichlet) eigenvalue problem in 1D to illustrate the idea and then use the tensor-product structure to generate the stiffness and mass matrices for multiple dimensional problems. To remove the outliers, we penalize the product of the higher-order derivatives from both the solution and test spaces at the domain boundary. Intuitively, we construct a better approximation by weakly imposing features of the exact solution. Effectively, we add terms to the variational formulation at the boundaries with minimal extra computational cost. We then generalize the idea to remove the outliers for the isogeometric analysis to the Neumann eigenvalue problem (for p ≥ 2). The boundary penalization does not change the test and solution spaces. In the limiting case when the penalty goes to infinity, we perform the dispersion analysis of C 2 cubic elements for Dirichlet eigenvalue problem and C 1 quadratic elements for Neumann eigenvalue problem. We obtain the analytical eigenpairs for the resulting matrix eigenvalue problems. Numerical experiments show optimal convergence rates for the eigenvalues and eigenfunctions of the discrete operator.