How to programmatically calculate real eigenvalues and (optionally) complex eigenvectors for big matrix I am trying to solve 1D timeindependent schrodinger equation
$$
-\frac{\hbar^{2}}{2 m} \frac{\partial^{2} \psi}{\partial x^{2}}+V(x) \psi=E \psi
$$
for periodic potential to simulate crystal lattice using computational mathematics methods:
$$\dfrac{\partial^2\psi}{\partial x^2}=\dfrac{\psi_{i+1}-2\psi_i+\psi_{i-1}}{\Delta x^2}$$
As I understand, I have to find eigenvalues for this matrix
$$
\left[\begin{array}{cccc}
\frac{1}{\Delta y^{2}}+L^{2} m V_{1} & -\frac{1}{2 \Delta y^{2}} & 0 & 0 \ldots \\
-\frac{1}{2 \Delta y^{2}} & \frac{1}{\Delta y^{2}}+L^{2} m V_{2} & -\frac{1}{2 \Delta y^{2}} & 0 \ldots \\
\ldots & \ldots & \ldots & -\frac{1}{2 \Delta y^{2}} \\
\ldots 0 & 0 & -\frac{1}{2 \Delta y^{2}} & \frac{1}{\Delta y^{2}}+L^{2} m V_{N-1}
\end{array}\right]\left[\begin{array}{cc}
\psi_1\\
\psi_2\\
\ldots\\
\psi_{N-1}\end{array}\right]=L^2mE\left[\begin{array}{cc}
\psi_1\\
\psi_2\\
\ldots\\
\psi_{N-1}\end{array}\right];\quad L^2=2/\hbar^2
$$
Perturbation theory for weak connection gives continuous spectrum, so probably I have to take a really big matrix to get characteristic polynomial of huge order. I am looking for methods or at least python library to solve such a big equation.
If it is important, I have nearly 5 GB of RAM and AMD® Athlon gold 3150u CPU, but probably I'll be able to use my university server
 A: Ok, I don't need to calculate complex eigenvectors because all eigenstates are independent and are real in some point in time. So with numpy for potential
$$U=\dfrac{U_0}{10(cos(i-5)\pi/5)^2+1};\quad U_0\approx 1\;eV$$
where $U_0$ I took nearly $1,6$ eV because experimentally silicium has 1,21 eV conduction band and where $i-$number if node (10 nodes per atom and 100 atoms).
I used this code to calculate eigenvectors
import numpy as np
import math
import scipy.linalg as la
from numpy.linalg import eigvalsh
def potential(i):
    return -1/(10*(math.cos((i-5)*math.pi/5)-1)*(math.cos((i-5)*math.pi/5)-1)+1)
def element(i,j):
    if i==j:
        return 0.1524+potential(i)
    if abs(i-j)==1:
        return -0.1524/2
    return 0
I = 1j
a=[]
for i in range(0,5000):
    b=[]
    for j in range(0,5000):
        b.append(element(i,j))
    a.append(b)
arr = np.array(a)
print('matrix created')
eigvals, eigvecs = la.eig(arr)
print('eigenvalues and eigenvectors calculated')
eigvals = eigvals.real
f = open('Energies.txt', 'w')
for element in eigvals:
    f.write(str(element)+'\n')
print("eigenvalues stored")
f = open('Psi_functions.txt', 'w')
for element in eigvecs:
    f.write(str(element)+'\n') 

If anybody need Psi functions for selfpurposes you can get them here.
It's much simplier to calculate eigenvalues so I made 200 nodes per atom:
$$U=\dfrac{U_0}{10(cos(i/20-5)\pi/5)^2+1};\quad U_0\approx 1\;eV$$
and used eigvalsh(arr) instead off la.eig(arr). So I got these energies:

