You're looking for a set of wavefunctions $\psi_n(x)$ that are solutions to $H\psi_n=E\psi_n$. The index $n$ does not refer to a discretized position. For numerical solutions, you will have to discretize $x$, for example with constant steps $h$ (you call them $\Delta x$); then the second derivative will be approximated as $$ {\partial^2\psi\over\partial x^2} \approx {\psi(x+h)+\psi(x-h)-2\psi(x) \over h^2}. $$ You can write $\psi_{nm}\equiv \psi_n(mh)$ and express the partial derivative above as a matrix multiplication. Make sure that $h$ is much smaller than the smallest expected wavelength of your wavefunctions. So your discretized Schrödinger equation becomes $$ {\hbar^2\over 2m} \sum_m \mathcal D_{jm}\psi_{nm} = V_j - E_n, $$$$ -{\hbar^2\over 2m} \sum_m \mathcal D_{jm}\psi_{nm} + V_j\psi_{nj} = E_n\psi_{nj}, $$ where $\mathcal D$ is the discretized differential operator. For example, if you have discretized your $x$ space into four possible $x$ values, $\{0,h,2h,3h\}$, then the differential operator is $$ \mathcal D = \frac{1}{h^2} \left( \begin{array}{cccc} -2 & 1 & 0 & 0 \\ 1 & -2 & 1 & 0 \\ 0 & 1 & -2 & 1 \\ 0 & 0 & 1 & -2 \end{array}\right). $$ In this representation, you have to deal carefully with the edges of the domain. You'd better choose your domain large enough that $\psi$ is very close to zero at those edges. The potential operator is a diagonal matrix, $$ V = \left( \begin{array}{cccc} V(0) & 0 & 0 & 0 \\ 0 & V(h) & 0 & 0 \\ 0 & 0 & V(2h) & 0 \\ 0 & 0 & 0 & V(3h) \end{array}\right) $$ And the total Hamiltonian is $$H_{jm}=-(h^2/2m)\mathcal D_{jm} + V_{jm},$$ of which you want to find eigenvectors.