The Schrödinger equation for a particle in a one-dimensional Gaussian potential well , given by , has never been solved analytically. This Demonstration derives an approximation for the first few bound-state energies, , using the linear variational method. The wavefunction is approximated by a linear combination . It is convenient to take the basis functions as the corresponding orthonormalized eigenfunction of the linear harmonic oscillator: , where is the Hermite polynomial and is a scaling constant to be determined variationally. After evaluating the matrix elements over the selected set of basis functions, Mathematica can calculate the eigenvalues in a single step, from which we select only those with negative values. For convenience, we set , so that all distances are expressed in bohrs (Bohr radii) and energy quantities in hartrees.