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.