The time-independent Schrodinger Equation is
$$ \frac{\hbar^2}{2m}\frac{d^2 \psi}{dx^2} + V(x)\psi = E\psi. $$
Now suppose we have a system with a given potential $V$ inside an infinite square well of length $L$. This gives the boundary conditions: $\psi(0) = \psi(L) = 0$.
We can transform the Schrodinger Equation in the following way
$$ \begin{align*} &\hbar = 1 \quad&&\text{Natural Units}\\ &y = \frac{x}{L}\quad&&\text{To switch to dimensionless units} \end{align*} $$
These transformations yields:
$$ -\frac{1}{2}\frac{d^2 \psi}{dy^2} + mL^2 V(y)\psi = mL^2 E\psi $$
where the boundary condition becomes: $\psi(y=0) = \psi(y=1) = 0$.
Discretising the second derivative using the Finite Difference Method gives:
$$ \frac{d^2 f}{dx^2} = \frac{f_{j+1}-2f_j+f_{j-1}}{\Delta x^2} $$
so the Schrodinger Equation becomes
$$ -\frac{1}{2}\frac{\psi_{j+1}-2\psi_j+\psi_{j-1}}{\Delta y^2} +mL^2V_j\psi_j = mL^2E\psi_j $$
which becomes
$$ -\frac{1}{2 \Delta y^2}\psi_{j+1} + \left(\frac{1}{\Delta y^2}+mL^2V_j \right)\psi_j - \frac{1}{2 \Delta y^2}\psi_{j-1} = mL^2E\psi_j. $$
This is essentially just a large system of equations spanning from $j=1$ to $j=N$. Therefore, we can represent this system using a tridiagonal matrix:
$$ \begin{bmatrix}\frac{1}{\Delta y^2}+mL^2V_1 & -\frac{1}{2 \Delta y^2} & 0 & 0...\\ -\frac{1}{2 \Delta y^2} & \frac{1}{\Delta y^2}+mL^2V_2 & -\frac{1}{2 \Delta y^2} & 0... \\ ...& ... & ... & -\frac{1}{2 \Delta y^2}\\...0 & 0 & -\frac{1}{2 \Delta y^2} & \frac{1}{\Delta y^2}+mL^2V_{N-1} \\ \end{bmatrix} \begin{bmatrix} \psi_1 \\ \psi_2 \\ ... \\ \psi_{N-1} \end{bmatrix} = mL^2 E \begin{bmatrix} \psi_1 \\ \psi_2 \\ ... \\ \psi_{N-1} \end{bmatrix} $$
This eigenvalue equation is the discrete version of Schrodinger’s Equation and can be solved numerically. We will now attempt it using python.
Lets plot the familiar quantum harmonic oscillator potential with