We present a new method for computing the eigenfunctions of a linear differential operator such that they satisfy a discrete orthogonality constraint, and can thereby be used to efficiently approximate physical measurement data. Classic eigenfunctions, such as those of Sturm-Liouville problems, are typically used for approximating continuous functions, but are cumbersome for approximating discrete data. We take a matrix-based variational approach to compute eigenfunctions based on physical problems which can be used in the same manner as the DCT or DFT. The approach uses sparse matrix methods, so it can be used to compute a small number of eigenfunctions. The method is verified on common eigenvalue problems, and the approximation of real-world measurement data of a bending beam by means of its computed mode-shapes.