  The conventional toy-model constructions of phase diagrams often use various
versions of the standard Hermitian Bose-Hubbard Hamiltonians $H$. These studies
were recently extended to cover several non-Hermitian PT-symmetric versions of
the model. A key technical merit of the extensions has been found in the
existence of Kato's phase-transition-related exceptional-point spectral
degeneracies of the $N$th order (EPN). In this setting, an ill-conditioned
nature of $H$ near its EPN singularity represented a decisive technical
obstacle. In our paper we show how the obstacle can be circumvented, leading to
a broader class of the eligible models with tridiagonal and symmetric complex
Hamiltonians $H^{(N)}(z)$ exhibiting the EPN degeneracies. The implementation
of the method (requiring the use of computer-assisted symbolic manipulations at
the larger integers $N$) is illustrated by sampling the real or complex spectra
of energies near as well as far from the EPN dynamical regime.
