I am trying to find a general method to discretize a finite crystal system.
How I have been discretizing systems so far (using Wannier functions):
When you have an infinite crystal, you may apply Bloch theorem to conclude your energy eigenfunctions are also eigenfunctions of crystal-translation operators, in which case they can be written as $\psi_{\lambda,\vec{k}}(\vec{r})=\exp\left(i\vec{k}\cdot\vec{r}\right)u_{\lambda,\vec{k}}\left(\vec{r}\right)$ where $u_{\lambda,\vec{k}}\left(\vec{r}\right)$ have the periodicity of the Hamiltonian: $u_{\lambda,\vec{k}}\left(\vec{r}\right)=u_{\lambda,\vec{k}}\left(\vec{r}+\vec{R}\right)$ where $\vec{R}$ are the lattice vectors, and $\psi_{\lambda,\vec{k}}(\vec{r})$ is periodic in $\vec{k}$: $\psi_{\lambda,\vec{k}}(\vec{r}) = \psi_{\lambda,\vec{k}+\vec{K}}(\vec{r})$ where $\vec{K}$ are the reciprocal lattice vectors to $\vec{R}$.
Due to the periodicity in $\vec{k}$, we may expand $\psi_{\lambda,\vec{k}}(\vec{r})$ in a Fourier series as $\psi_{\lambda,\vec{k}}(\vec{r}) = \sum_{\vec{R}} \exp\left(i\vec{k}\cdot\vec{R}\right)\phi_{\lambda,\vec{R}}\left(\vec{r}\right)$. The $\phi_{\lambda,\vec{R}}\left(\vec{r}\right)$ functions are called "Wannier" functions and they are orthogonal (between different bands $\lambda$ and different values of $\vec{R}$). Using this decomposition we are able to write a general wave-function $\Psi\left(x\right)$ as \begin{align} \Psi\left(x\right)&=\sum_{\lambda,\vec{k}}\alpha_{\lambda,\vec{k}}\psi_{\lambda,\vec{k}}(\vec{r})\\&=\sum_{\lambda,\vec{k}}\alpha_{\lambda,\vec{k}}\sum_{\vec{R}}\exp\left(i\vec{k}\cdot\vec{R}\right)\phi_{\lambda,\,\vec{R}}\left(\vec{r}\right)\\&=\sum_{\vec{R}}\sum_{\lambda}\underbrace{\left(\sum_{\vec{k}}\alpha_{\lambda,\vec{k}}\exp\left(i\vec{k}\cdot\vec{R}\right)\right)}_{\beta_{\lambda,\,\vec{R}}}\phi_{\lambda,\,\vec{R}}\left(\vec{r}\right)\\&=\sum_{\vec{R}}\sum_{\lambda}\beta_{\lambda,\,\vec{R}}\phi_{\lambda,\,\vec{R}}\left(\vec{r}\right) \end{align}
As a result, we can work with $\beta_{\lambda,\,\vec{R}}$ as wave-functions that depend on a discrete parameter ($\lambda, \vec{R}$) instead of $\Psi\left(x\right)$ and write everywhere (infinite) matrix equations instead of differential equations for the Schroedinger equation (by writing matrix-elements with the Wannier functions).
For example, instead of writing $$H\Psi\left(x\right)=E\Psi\left(x\right)$$ (a differential equation) we can write $$\sum_{\lambda',\vec{R}'}H_{\lambda,\vec{R},\lambda',\vec{R}'}\beta_{\lambda',\,\vec{R}'} = E \beta_{\lambda,\,\vec{R}} \tag{2}$$ (an (infinite) matrix equation), where $H_{\lambda,\vec{R},\lambda',\vec{R}'} \equiv \int_{\mathbb{R}^3}\phi_{\lambda,\vec{R}}\left(\vec{r}\right)^\ast H\phi_{\lambda,\vec{R}}\left(\vec{r}\right)d^3\vec{r}$.
My question Is
How do I convert the Schroedinger equation to an equation of the form $\left(2\right)$ in situations when I cannot apply the Bloch theorem (and so may not build Wannier functions)? Such instances could be: (a) the existence of a perturbing potential which does not have the periodicity of the potential or (b) the finite nature of the crystal where we do not take to have periodic boundary conditions but actual edge boundary conditions (such as here).
I am looking for a general way to do this process.
I have read some articles that lead me to pessimistic conclusions that this might be complicated, for example Wilkinson (1998).
