  We present higher order polynomial algebras which are the dynamical symmetry
algebras of a wide class of multi-mode boson systems in non-linear optics. We
construct their unitary representations and the corresponding single-variable
differential operator realizations. We then use the results to obtain exact
(Bethe ansatz) solutions to the multi-mode boson systems, including the BEC
models as special cases.
