We introduce a simple, general, and convergent scheme to compute generalized eigenfunctions of self-adjoint operators with continuous spectra on rigged Hilbert spaces. Our approach does not require prior knowledge about the eigenfunctions, such as asymptotics or other analytic properties. Instead, we carefully sample the range of the resolvent operator to construct smooth and accurate wave packet approximations to generalized eigenfunctions. We prove high-order convergence in key topologies, including weak-star convergence for distributional eigenfunctions, uniform convergence on compact sets for locally smooth generalized eigenfunctions, and convergence in seminorms for separable Frechet spaces, covering the majority of physical scenarios. The method's performance is illustrated with applications to both differential and integral operators, culminating in the computation of spectral measures and generalized eigenfunctions for an operator associated with Poincare's internal waves problem. These computations corroborate experimental results and highlight the method's utility for a broad range of spectral problems in physics.