  We present a simple construction for a tridiagonal matrix $T$ that commutes
with the hopping matrix for the entanglement Hamiltonian ${\cal H}$ of open
finite free-Fermion chains associated with families of discrete orthogonal
polynomials. It is based on the notion of algebraic Heun operator attached to
bispectral problems, and the parallel between entanglement studies and the
theory of time and band limiting. As examples, we consider Fermionic chains
related to the Chebychev, Krawtchouk and dual Hahn polynomials. For the former
case, which corresponds to a homogeneous chain, the outcome of our construction
coincides with a recent result of Eisler and Peschel; the latter cases yield
commuting operators for particular inhomogeneous chains. Since $T$ is
tridiagonal and non-degenerate, it can be readily diagonalized numerically,
which in turn can be used to calculate the spectrum of ${\cal H}$, and
therefore the entanglement entropy.
