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$; 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, $$
where $\mathcal D$ is the discretized differential operator.