In this seminar, we delve into the diagonalization of the so-called Heun-Askey-Wilson (HAW) operator. To achieve this, we employ the modified algebraic Bethe ansatz technique along with the theory of Leonard pairs. The (dual) eigenstates of the HAW operator are constructed as (dual) Bethe states, where the (dual) Bethe roots satisfy Bethe ansatz equations of homogeneous or inhomogeneous type. A precise correspondence between (dual) eigenvectors of elements of a Leonard pair of q-Racah type and (dual) Bethe states can be established. This leads to various expressions of the q-Racah polynomials in terms of ratios of scalar products of Bethe states. We then discuss the possibility of expressing the q-Racah polynomials in terms of Slavnov-Gaudin-Korepin determinants.