In this thesis, p. 27, the following Euler equation is given, see (4.9): $$ u_{i+1,j}=\frac{h_t}{2h_x^2}\left[\left(\frac{2h_x^2}{h_t}+2a_1h_x^2-4a_2\right)u_{i,j}+(2a_2-a_3h_x)u_{i,j-1}+(2a_2+a_3h_x)u_{i,j+1}+2b_1h_x^2\right] $$
Here, $i$ represents the time and $j$ represents the (1d-) space. Moreover, $h_t$ and $h_x$ are the step sizes with respect to time and space, respectively. $a_1,a_2$ and $b_1$ are constants.
Then, on both sides, the z-tranformation (with respect to time) is applied, where $$ \mathcal{Z}(u_{i+1,j})=:U_j,~\mathcal{Z}(u_{i,j})=z^{-1}U_j,~\mathcal{Z}(u_{i,j-1})=z^{-1}U_{j-1},~\mathcal{Z}(u_{i,j+1})=U_{j+1}. $$
Due to the linked thesis, this gives (4.10) $$ U_j=\frac{h_t}{2h_x^2}\left[\left(\frac{2h_x^2}{h_t}+2a_1h_x^2-4a_2\right)z^{-1}U_j+(2a_2-a_3h_x)z^{-1}U_{j-1}+(2a_2+a_3h_x)z^{-1}U_{j+1}+\color{blue}{2b_1h_x^2}\right] $$
I am really wondering about the blue summand! Why isn't it $$ \mathcal{Z}(2b_1h_x^2)=2b_1h_x^2\mathcal{Z}(1)? $$
Is there some reason for that or is it just a mistake? I would really prefer that there is some reason for it. :-)
In particular this seems relevant since the author is interested in poles and zeros $z$ of the z-transform with $\lvert z\rvert <1$ for stability reasons; but for $\lvert z\rvert <1$, we have that $\mathcal{Z}(1)$ diverges. Moreover, the determination of the zeros and poles does not work as done in the thesis in case the summand is $2b_1h_x^2\mathcal{Z}(1)$.
For the equation that you write, the Z-transform does indeed include the $\mathcal Z(1)$ term.
I think there should be a term with the initial condition $u_{0,j} $, because $\mathcal Z(u_{i,j}) = z^{-1}U_j + u_{0,j}$, but I digress.
Here is an idea to get rid of having to deal with $b_1$, at least in case $a_1 \ne 0$, because in that case, the constant function $v_{i,j} = -\frac{b_1}{a_1}$ satisfies the non-homogeneous equation, and so, $u-v$ satisfies the homogeneous equation (without $b_1$).