The spectrum of the Laplacian operator can only be computed analytically for a few simple geometries. When direct computation is not possible, numerical methods can be useful for finding approximations of the spectrum. In this paper, we present two major techniques for approximating the spectrum of the Laplacian operator via discretization, namely the finite difference and finite element methods. We show how to construct the matrices whose eigenvalues approximate the desired spectrum, and we discuss the convergence of these methods as the grid is refined. We also introduce a few numerical methods that are commonly used for solving the associated matrix eigenvalue problems. Finally, we illustrate how the above techniques can be used to handle higher order operators, such as the biharmonic operator.