Simple, combine both real- and $\mathbf{k}$-space pictures! The basic idea is to split up your $n$-dimensional system into multiple $(n-1)$-dimensional systems. For example, say you have a 2D square lattice and you define your edges along the $x$-direction. Then you need to break the 2D lattice into 1D lattices pointing in the $x$-direction. In other words, you need to break translational symmetry in the $y$-direction. For the sake of (analytical) simplicity, consider the model discussed in:
Markus König, Hartmut Buhmann, Laurens W. Molenkamp, Taylor Hughes, Chao-Xing Liu, Xiao-Liang Qi, and Shou-Cheng Zhang. “The Quantum Spin Hall Effect: Theory and Experiment.” Journal of the Physical Society of Japan 77, no. 3 (2008): 031007. (arXiv)
In Eq. (10) they have a $\mathbf{k}$-space model, also known as the Bernevig-Hughes-Zhang (BHZ) model, of the entire 2D system $${\cal H}=\sum_{\mathbf{k}}\left(A\sin(k_{x})\Gamma^{1}+A\sin(k_{y})\Gamma^{2}+{\cal M}(\mathbf{k})\Gamma^{5}\right)c_{\mathbf{k}}^{\dagger}c_{\mathbf{k}}$$ where $M(\mathbf{k})=M-2B\left(2-\cos(k_{x})-\cos(k_{y})\right)$ and the lattice constant has been set to 1. The next step, as Eq. (11) indicates, is to Fourier transform back to real-space in only the $y$-direction but leave the $x$-direction unchanged. That's what I meant when I said we “combine real- and $\mathbf{k}$-space pictures.” In other words, we are breaking translational symmetry only in the $y$-direction. This is done by plugging Eq. (11), which repeated here for convenience, into the above equation $$c_{\mathbf{k}}=\frac{1}{L}\sum_{j}e^{ik_{y}j}c_{k_{y},j}$$ where $j$ is the lattice (or 1D chain) coordinate in the $y$-direction and $y=0,1,2,\dots L$. In this example $L$ will not matter as much; but it will when you compute the dispersion numerically. A brute force plug-and-play gives
$$\mathcal{H} = \frac{1}{L^{2}}\sum_{k_{x}k_{y}}\left[A\sin(k_{x})\Gamma^{1}+\frac{A}{2i}\left(e^{ik_{y}}-e^{-ik_{y}}\right)\Gamma^{2}\right. \\ \left.+\left(M-2B\left(2-\cos(k_{x})-\frac{1}{2}\left(e^{ik_{y}}+e^{-ik_{y}}\right)\right)\right)\Gamma^{5}\right]\sum_{\ell}e^{-ik_{y}\ell}c_{k_{x},\ell}^{\dagger}\sum_{j}e^{ik_{y}j}c_{k_{x},j}$$
$$\mathcal{H} = \frac{1}{L}\sum_{k_{x}}\sum_{j}\left[A\sin(k_{x})\Gamma^{1}+\left(M-4B+2B\cos(k_{x})\right)\Gamma^{5}\right]c_{k_{x},j}^{\dagger}c_{k_{x},j} \\
+\frac{1}{L}\sum_{k_{x}}\sum_{j}\left(-\frac{iA}{2}\Gamma^{2}+B\Gamma^{5}\right)c_{k_{x},j+1}^{\dagger}c_{k_{x},j}$$
$$\mathcal{H} = \frac{1}{L}\sum_{k_{x}}\sum_{j}\mathcal{M}(k_{x})c_{k_{x},j}^{\dagger}c_{k_{x},j}+\frac{1}{L}\sum_{k_{x}}\sum_{j}\mathcal{T}^{\dagger}c_{k_{x},j+1}^{\dagger}c_{k_{x},j}+\frac{1}{L}\sum_{k_{x}}\sum_{j}\mathcal{T}c_{k_{x},j-1}^{\dagger}c_{k_{x},j}$$
where $\mathcal{M}(k_{x})=A\sin(k_{x})\Gamma^{1}+\left(M-4B+2B\cos(k_{x})\right)\Gamma^{5}$ and $\mathcal{T}=(iA/2)\Gamma^{2}+B\Gamma^{5}$ and we have made use of the delta function identity of the type $$\frac{1}{L}\sum_{k_{y}}e^{ik_{y}(j-\ell\pm1)}=\delta_{j-\ell\pm1}$$ several times. With the ansatz in Eq. (15), i.e. $\psi_{\alpha}(j)=\lambda^{j}\phi_{\alpha}$, an analytic solution of the edge states can be obtained (see Eq. (22)). The solution of the eigenvalue equation using this ansatz has been elegantly discussed in section 2.2 of:
G. Tkachov, and E. M. Hankiewicz. “Spin-helical transport in normal and superconducting topological insulators.” physica status solidi (b) 250, no. 2 (2013): 215. (arXiv)
and I will not repeat it here. In the case of graphene, as discussed by Kane and Mele, we are not so fortunate. In that case, we need to diagonalize the above Hamiltonian numerically by choosing $L$ = 50-100. The main criterion in determining $L$ is making sure that the edge state wave function overlap at opposite boundaries ($y$=0 and $y=L$) is negligible. My guess is that you just figure it out by trial and error.
Another main difference between BHZ and the Kane-Mele model is that in the Kane-Mele model we have the added complexity of determining whether we have a zig-zag or armchair boundary. Depending on what choice we make, we need to define the 1D systems accordingly; they obviously won't be straight lines, as in BHZ, and will depend on whether you break translational symmetry in the $x$- or $y$-direction.
Hope that helped.
PS: I know I have skipped a bunch of steps in the above algebraic manipulations and referred the rest of the solution to the above paper. In case you're interested I could upload a PDF document containing all the steps.
