This answer, I hope, complements those earlier ones and give some little explanation as to how the basis set calculation works.
As pointed out already it is not necessary to use a basis set to solve the Schroedinger equation because for simple systems, such as a harmonic or Morse potential, analytical i.e. algebraic solutions are available. Numerical methods, such as the Numerov and Shooting methods, can be used for more complicated potentials but soon become intractable.
Generating a solution using a basis set is often computationally the best method, and sometimes the only method that can be used. This method involves solving the Schroedinger equation without direct integration.
The idea is to choose a known set of functions which when added together in the correct proportions will generate the wavefunction for the target problem. This is the same idea as using a fourier series to generate a target function; the more terms that are included in the series the better the fit to the target function becomes. An example of a fourier series is to form a polynomial using many sine/cosine waves of different amplitude and frequency.
In forming a basis set we choose a set of orthogonal functions, and numerically solve the matrix Schroedinger equation using as many of these functions as needed for the accuracy required. Normalised orthogonal functions have the property that $\int \phi^*_n\phi_mdx=\delta_{nm}$ where $\delta_{nm}$ is the Knonecker delta, if $n=m$ then $\delta_{nm}=1$ otherwise it is zero.
As an example we may use the sine-wave solutions to the particle in a box (PIB) problem to solve a harmonic well potential with a central 'hill' or barrier. If the PIB equation is written as $H\phi = E\phi$, the eigenvalues (energies) are
$\epsilon_n=\int\phi^*_n~H~\phi_ndx$.
The new problem we want to solve has the equation
$$ H^{’} \psi = E^{’} \psi \tag {1} $$
which has a different Hamiltonian to reflect the different potential energy. The new wavefunction for the $n^{th}$ level is assumed to have the form
$$\psi_n=c_{1,n}\phi_1 + c_{2,n}\phi_2 \cdots c_{N,n}\phi_N$$
where the coefficients c are to be found. (The calculation follows that in J.Chem. Educ. 2012, v89, p1152). Perhaps twenty or a hundred terms in the basis set will be needed to obtain an accurate answer both for the eigenvalues (energies) and eigenvectors (wavefunctions). If $N \rightarrow \infty$ then the energy and wavefunctions would be exact.
The PIB basis set used here is a good starting point in our example as we know that the wavefunctions for a harmonic oscillator oscillate somewhat in a sinusoidal manner and the barrier harmonic well should be similar. If we used exponential functions many, many more would be needed to achieve the same result as these just do not have the same shape as the wavefunction. However, exponential functions would be good to use in describing an atoms wavefunctions as they have a shape similar to atomic orbital radial functions.
The energy for level m with the new Hamiltonian $H^{’}$ is found by left multiplying eqn. 1 by $\phi^*$
$$\int \phi^*_mH^{’} \psi_ndx =E^{’}_n\int \phi^*_m\psi_ndx$$
and substituting for $\psi$ gives
$$\int \phi^*_mH^{’}(c_{1,n}\phi_1 + c_{2,n}\phi_2 \cdots)dx=E^{’}_n\int \phi^*_m(c_{1,n}\phi_1 + c_{2,n}\phi_2 \cdots)dx$$
On expanding this out and using the orthonormality the equation reduces to a matrix equation
$$ \bf HC=E^{’}C$$
where $\bf H$ is a square matrix of elements
$$H_{m,n} = \int \phi^*_mH^{’}\phi_ndx$$
and $\bf C$ is a column matrix of coefficients. The matrix equation is $\bf(H-E^{’}I)C=0$ where $\bf I$ is the unit diagonal matrix and $\bf 0$ a null matrix containing only zeros. Solving this equation involves evaluating many simple integrals $H_{m,n}$ in the secular determinant and then diagonalising to find the characteristic equation which is a polynomial in $E^{’}$ and hence the energies.