  In this work we present a general construction of integrable models for boson
tunneling in multi-well systems. We show how the models may be derived through
the Quantum Inverse Scattering Method and solved by algebraic Bethe ansatz
means. From the transfer matrix we find only two conserved operators. However,
we construct additional conserved operators through a different method. As a
consequence the models admit multiple pseudovacua, each associated to a set of
Bethe ansatz equations. We show that all sets of Bethe ansatz equations are
needed to obtain a complete set of eigenstates.
