  We propose a robust and efficient algorithm for computing bound states of
infinite tight-binding systems that are made up of a finite scattering region
connected to semi-infinite leads. Our method uses wave matching in close
analogy to the approaches used to obtain propagating states and scattering
matrices. We show that our algorithm is robust in presence of slowly decaying
bound states where a diagonalization of a finite system would fail. It also
allows to calculate the bound states that can be present in the middle of a
continuous spectrum. We apply our technique to quantum billiards and the
following topological materials: Majorana states in 1D superconducting
nanowires, edge states in the 2D quantum spin Hall phase, and Fermi arcs in 3D
Weyl semimetals.
