Once I have the eigenvalues and the eigenvectors, how do I find the eigenfunctions? I am using Mathematica to construct a matrix for the Hamiltonian of some system. I have built this matrix already, and I have found the eigenvalues and the eigenvectors, I am uncertain if what I did next is correct: I took the normalized eigenvectors, placed them in matrix form, and did matrix multiplication with the basis set of solutions.
Let me try to be more precise since I am not sure I am using the right language when mentioning the basis solutions. In the problem we are using the set of solutions of the particle in a box model as our basis. I can increase the number of basis elements in the calculation of the matrix of the Hamiltonian (which amounts to doing $<\psi_n|H \psi_k>$ over a specified range of $n$ and $k$) in order for some of my smallest eigenvalues to begin to converge. Once I have this $H$ matrix built, and that I see that my eigenvalues are converging to some degree, I take the eigenvectors of the $H$ matrix, format them to be in matrix form, and multiply them by the set of basis solutions.
I hope that makes things clearer.
 A: If $\bf{v}$ is an eigenvector of the matrix $\bf{H}$ (where the ith row and jth column of $\bf{H}$ is $<\psi_i|H|\psi_j>$) with eigenvalue $\lambda$, i.e.
$$\bf{H} \cdot \bf{v} = \lambda \cdot \bf{v}$$
the function (which is the one you are looking for)
$$ \varphi = \bf{v} \cdot \bf{\psi} = \sum_j \bf{v}_j \cdot \psi_j$$
is an 'eigenfunction' (solution of the Schrödinger equation) of the Hamiltonian corresponding to $\bf{H}$ because for all $i$:
$$ <\psi_i|H|\varphi> = \sum_j \bf{v}_j \cdot <\psi_i|H|\psi_j> = \sum_j \bf{v}_j \cdot \bf{H}_{ij} = (\bf{H} \cdot \bf{v})_i = (\lambda \cdot \bf{v})_i$$
Assuming your set of basis functions is orthonormal, i.e.
$<\psi_i|\psi_j>= \delta_{ij}$ one can rewrite the above expression as:
$$<\psi_i|H|\varphi> = \sum_j \delta_{ij} \cdot \lambda \cdot \bf{v}_j = \lambda \cdot \sum_j \bf{v}_j <\psi_i|\psi_j> = \lambda \cdot <\psi_i|\sum_j \psi_j> = \lambda \cdot <\psi_i|\varphi>$$
because this holds for all $i$
$$H|\varphi> = \lambda \cdot |\varphi>$$
holds.

You say that you put the eigenvector $\bf{v}$ in Matrix form and then multiply it with the vector of basis functions to obtain the function $\varphi$. In fact it should be more like a 'dot product' but if you put the numbers of the eigenvector onto the diagonal (and leave zeros off the diagonal), that should be equivalent.
