Discretization of the Schrödinger equation
$$\require{color}$$ The Hamiltonian matrix in the continous position state basis is given by:
$$H(x,x')=\int_{-\infty}^{+\infty}\text{d}p\frac{p^{2}}{2m}\langle x|p\rangle\langle p|x'\rangle+\langle x|V(\hat{x})|x'\rangle$$
From the quantum mechanics we know that the scalar product of the position and momentum eigenstates is given by:
$$\langle x|p\rangle=\frac{1}{\sqrt{2\pi\hbar}}\exp\left(\frac{\text{i}px}{\hbar}\right),$$
being nothing else but the position representation of a momentum eigenstate. By taking the complex conjugate one can also replace the second scalar product by:
$$\langle p|x'\rangle=\frac{1}{\sqrt{2\pi\hbar}}\exp\left(-\frac{\text{i}px'}{\hbar}\right)$$
Thus, the matrix element can be rewritten as:
$$H(x,x')=\frac{1}{2\pi\hbar}\int_{-\infty}^{+\infty}\text{d}p\frac{p^{2}}{2m}\exp\left(\frac{\text{i}p(x-x')}{\hbar}\right)+\langle x|V(\hat{x})|x'\rangle$$
Already here, we recognize a continous Fourier transform in the kinetic part of the hamiltonian expression!
We can now discretize this by introducing a grid in the position and the momentum space :
$$x_{i}=x_{min}+i\Delta x\qquad i=0,\cdots,N-1$$
$$p_{k}=p_{min}+k\Delta p\qquad k=0,\cdots N-1$$
In this way, the continous matrix \( H(x,x') \) is represented by a discrete matrix:
$$H(x,x')\rightarrow h_{ij}=\frac{1}{2\pi\hbar}\sum_{k=0}^{N-1}\Delta p\frac{p_{k}^{2}}{2m}\exp\left(\frac{\text{i}p_{k}(i-j)\Delta x}{\hbar}\right)+\frac{1}{\Delta x}V(x_{i})\delta_{ij}$$
Here, we have utilized the fact that:
$$\langle x|V(\hat{x})|x'\rangle = V(x)\langle x|x'\rangle = V(x)\delta(x-x')$$
which can be transformed into
$$\langle x|V(\hat{x})|x'\rangle\rightarrow\langle x_{i}|V(\hat{x})|x_{j}\rangle=V(x_{i})\delta\left[(i-j)\Delta x\right]=\frac{1}{\Delta x}V(x_{i})\delta_{ij}$$
by utilizing the scaling property of the Dirac's delta function:
$$\delta (\alpha x)={\frac {\delta (x)}{|\alpha |}}.$$
By utilizing the positon space grid, the continous version of the eigenvalue problem can be turned into:
$$\begin{array}{c} \int_{-\infty}^{+\infty}\text{d}x'H(x,x')\psi(x')=E\psi(x)\\ \downarrow \end{array}$$
$$\sum_{j=1}^{N-1}\left[\sum_{k=0}^{N-1}\frac{\Delta x\Delta p}{2\pi\hbar}\frac{p_{k}^{2}}{2m}\exp\left(\frac{\text{i}p_{k}(i-j)\Delta x}{\hbar}\right)+V(x_{i})\delta_{ij}\right]\psi_{j}=E\psi_{i}\quad i=0,\cdots,N-1$$
where \( \psi_{i} :=\psi(x_{i}) \) and \( \psi_{j} :=\psi(x_{j}) \) has been introduced as a shorthand notation.
This now represents an usual matrix eigenvalue problem:
$$\sum_{j=0}^{N-1}H_{ij}\psi_{j}=E\psi_{i}\quad i=0,\cdots,N-1,$$
with hamiltonian matrix elements defined by:
$$H_{ij}=\sum_{k=0}^{N-1}\frac{\Delta x\Delta p}{2\pi\hbar}\frac{p_{k}^{2}}{2m}\exp\left(\frac{\text{i}p_{k}(i-j)\Delta x}{\hbar}\right)+V(x_{i})\delta_{ij}$$
In order to finish the discretization, we will now transform the hamiltonian matrix elements in a shape that will allow us to calculate them using standard Fast-Fourier-Transform (FFT) routine, which is available in the numpy.fft package.
We recall that:
$$ p_{k}=p_{min}+k\Delta p\quad k=0,\cdots,N-1 $$
This can be inserted into the expression for the matrix element: $$H_{lj}=\sum_{k=0}^{N-1}\frac{\Delta x\Delta p}{2\pi\hbar}\frac{\left[p_{min}+k\Delta p\right]^{2}}{2m}\exp\left(\frac{\text{i}(p_{min}+k\Delta p)(i-j)\Delta x}{\hbar}\right)+V(x_{i})\delta_{ij}$$
and can be simplified to:
$$H_{ij}=\exp\left(\frac{\text{i}(p_{min}(i-j)\Delta x}{\hbar}\right)\sum_{k=0}^{N-1}\frac{\Delta x\Delta p}{2\pi\hbar}\frac{\left[p_{min}+k\Delta p\right]^{2}}{2m}\exp\left(\frac{\text{i}k(i-j)\Delta x\Delta p}{\hbar}\right)+V(x_{i})\delta_{ij}$$
Since the integration range in the momentum space has to be symmetric, the lowest and the highest value of the discretized momentum need to fullfil the following relation:
$$p_{0}=p_{min}$$
and
$$p_{N-1}=-p_{min}-\Delta p=p_{min}+(N-1)\Delta p$$
This can be exploited to calculate $p_{min}$ and eliminate it from the final expression for the discretized momentum. From the above boundary values we obtain:
$$p_{min}=-\left(\frac{N}{2}\right)\Delta p $$
which can be inserted into the matrix element leading to:
$$ H_{ij}=\exp\left(\frac{-\text{i}N(i-j)\Delta x\Delta p}{2\hbar}\right)\sum_{k=0}^{N-1}\frac{\Delta x\Delta p}{2\pi\hbar}\frac{\left[p_{min}+k\Delta p\right]^{2}}{2m}\fcolorbox{red}{#f2d5d3}{$\exp\left(\frac{\text{i}k(i-j)\Delta x\Delta p}{\hbar}\right)$}+V(x_{i})\delta_{ij} $$
In the next step, we transform the sum over $k$ into a discrete Fourier transform. By comparing the red highlighted term in the above equation with the standard definition of the DFT:
$$ F_{k}=\sum_{j=0}^{N-1}f_{j}\exp\left[-\frac{2\pi \text{i} kj}{N}\right]\quad k=1,\cdots N-1. $$
we see, that by imposing the condition:
$$\frac{\Delta x\Delta p}{\hbar}=\frac{2\pi}{N}$$
the sum in the above matrix element will start to resemble a discrete Fourier transfrom:
$$H_{ij}=\frac{1}{N}\exp\left(-\text{i}\pi(i-j)\right)\sum_{k=0}^{N-1}\frac{\left[\frac{2\pi\hbar}{N\Delta x}\left(k-\frac{N}{2}\right)\right]^{2}}{2m}\exp\left(\frac{2\pi \text{i}ki}{N}\right)\exp\left(-\frac{2\pi \text{i}kj}{N}\right)+V(x_{i})\delta_{ij}$$
This can be seen more explicitly by defining a vector with the components:
$$f_{ik}=\frac{\left[\frac{2\pi\hbar}{N\Delta x}\left(k-\frac{N}{2}\right)\right]^{2}}{2m}\exp\left(\frac{2\pi \text{i}ki}{N}\right)$$
so that the hamiltonian matrix becomes:
$$H_{ij}=\exp\left(-\text{i}\pi(i-j)\right)\frac{1}{N}\sum_{k=0}^{N-1}f_{ik}\exp\left(-\frac{2\pi \text{i}kj}{N}\right)+V(x_{i})\delta_{ij}$$
or
$$H_{ij}=(-1)^{i-j}\frac{1}{N}\sum_{k=0}^{N-1}f_{ik}\exp\left(-\frac{2\pi \text{i}kj}{N}\right)+V(x_{i})\delta_{ij}$$