Given that $dr=(\eta-\gamma r)dt+\sqrt{\alpha r+\beta}dW$
Let $Z(r,t)=e^{A(t;T)-rB(t;T)}$,
\begin{matrix} \frac{dA}{dt}=\eta B-\frac{1}{2}\beta {{B}^{2}} \\ \frac{dB}{dt}=\frac{1}{2}\alpha {{B}^{2}}+\gamma B-1 \\ \end{matrix}
How can we solve for $B$ in this case?
$B$ is a Riccati Equation. The answer for B is given below:
$B(t;T)=\frac{2(e^{\Psi_1(T-t)}-1)}{(\gamma+\Psi_1)(e^{\Psi_1(T-t)}-1)+2\Psi_1}$,
where $\Psi_1=\sqrt{\gamma^2+2\alpha}$
As you noted, this is a Riccati type ODE and it can thus be simplified using the standard transformations for this class - see e.g. Wikipedia. We start by defining
\begin{equation} C(t, T) = \frac{1}{2} \alpha B(t, T) \qquad \Rightarrow \qquad C_t(t, T) = \frac{1}{2} \alpha B_t(t, t) \end{equation}
and get
\begin{eqnarray} C_t(t, T) & = & C^2(t, T) + \gamma C(t, T) - \frac{1}{2} \alpha. \end{eqnarray}
We now set
\begin{eqnarray} C(t, T) = -\frac{D_t(t, T)}{D(t, T)} \qquad \Rightarrow \qquad C_t(t, T) & = & -\frac{D_{tt}(t, T)}{D(t, T)} + \frac{D_t^2(t, T)}{D^2(t, T)}\\ & = & \frac{-D_{tt}(t, T)}{D(t, T)} + C^2(t, T) \end{eqnarray}
and get
\begin{eqnarray} \frac{D_{tt}(t, T)}{D(t, T)} & = & C^2(t, T) - C_t(t, T)\\ & = & -\gamma C(t, T) + \frac{1}{2} \alpha\\ & = & \gamma \frac{D_t(t, T)}{D(t, T)} + \frac{1}{2} \alpha \end{eqnarray}
or
\begin{equation} D_{tt}(t, T) = \gamma D_t(t, T) + \frac{1}{2} \alpha D(t, T). \end{equation}
This is a homogeneous second order linear ODE with constant coefficients and can be solved using standard methods. We note that $T$ has been fixed and make another substitution by defining $\tau = T - t$ such that $E(\tau) = D(t, T)$. We get
\begin{equation} E_{\tau \tau}(\tau) + \gamma E_\tau(\tau) - \frac{1}{2} \alpha E(\tau) = 0. \end{equation}
The characteristic equation is
\begin{equation} r^2 + \gamma r - \frac{1}{2} \alpha = 0 \end{equation}
with solutions
\begin{equation} r_{1, 2} = -\frac{1}{2} \gamma \pm\frac{1}{2} \sqrt{\gamma^2 + 2 \alpha} := \beta \pm \lambda. \end{equation}
Note that $\lambda = \frac{1}{2} \psi_1$ in your notation. We thus have the general solution
\begin{equation} E(\tau) = c_1 e^{(\beta + \lambda) \tau} + c_2 e^{(\beta - \lambda) \tau} \end{equation}
with
\begin{equation} E_\tau(\tau) = (\beta + \lambda) e^{(\beta + \lambda) \tau} + (\beta - \lambda) c_2 e^{(\beta - \lambda) \tau} \end{equation}
and for some constants $c_1$ and $c_2$ that have to be determined. We obtain the solution to the Riccati ODE by substituting back
\begin{equation} B(t, T) = \frac{2 C(t, T)}{\alpha} = -\frac{2 D_t(t, T)}{\alpha D(t, T)} = \frac{2 E_\tau(\tau)}{\alpha E(\tau)}. \end{equation}
Applying the terminal condition yields
\begin{eqnarray} B(T, T) = 0 \qquad & \Leftrightarrow & \qquad E_\tau(0) = 0\\ & \Leftrightarrow & \qquad c_1 = -c_2 \frac{\beta - \lambda}{\beta + \lambda}. \end{eqnarray}
Thus
\begin{equation} E(\tau) = \frac{c_2}{\beta + \lambda} e^{\beta \tau} \left( (\beta + \lambda) e^{-\lambda \tau} - (\beta - \lambda) e^{\lambda \tau} \right) \end{equation}
and
\begin{equation} E_\tau(\tau) = (\beta - \lambda) c_2 e^{\beta \tau} \left( e^{-\lambda \tau} - e^{\lambda \tau} \right). \end{equation}
We finally get
\begin{eqnarray} B(t, T) & = & \frac{2 (\beta^2 - \lambda^2) \left( e^{-\lambda \tau} - e^{\lambda \tau} \right)}{\alpha \left( (\beta + \lambda) e^{-\lambda \tau} - (\beta - \lambda) e^{\lambda \tau} \right)}\\ & = & \frac{\left( e^{2 \lambda \tau} - 1 \right)}{(\beta + \lambda) - (\beta - \lambda) e^{2 \lambda \tau}}\\ & = & \frac{2 \left( e^{\psi_1 (T - t)} - 1 \right)}{\left( \gamma + \psi_1 \right) \left( e^{\psi_1 (T - t)} - 1 \right) + 2 \psi_1}. \end{eqnarray}
