An operator of Heun-Askey-Wilson type is diagonalized within the framework of the algebraic Bethe ansatz using the theory of Leonard pairs. For different specializations and the generic case, the corresponding eigenstates are constructed in the form of Bethe states, whose Bethe roots satisfy Bethe ansatz equations of homogeneous or inhomogenous type. For each set of Bethe equations, an alternative presentation is given in terms of 'symmetrized' Bethe roots. Also, two families of on-shell Bethe states are shown to generate two explicit bases on which a Leonard pair acts in a tridiagonal fashion. In a second part, the (in)homogeneous Baxter T-Q relations are derived. Certain realizations of the Heun-Askey-Wilson operator as second q-difference operators are introduced. Acting on the Q-polynomials, they produce the T-Q relations. For a special case, the Q-polynomial is identified with the Askey-Wilson polynomial, which allows one to obtain the solution of the associated Bethe ansatz equations. The analysis presented can be viewed as a toy model for studying integrable models generated from the Askey-Wilson algebra and its generalizations. As examples, the q-analog of the quantum Euler top and various types of three-sites Heisenberg spin chains in a magnetic field with inhomogeneous couplings, three-body and boundary interactions are solved. Numerical examples are given. The results also apply to the time-band limiting problem in signal processing.MSc: 81R50; 81R10; 81U15.