I'm hoping to find some guidance on solving the following non-linear SDE in distribution. Meaning the distribution of $x_{t + \Delta t} | x_t$ where $x_t = [y_t, z_t]^\prime$.
$ \left[ \begin{array}{c} dy_t \\ dz_t \\ \end{array}\right] \sim \left[ \begin{array}{c} \theta_y \left( \alpha - y_t \right) + \zeta_0 z_t + \zeta_1 z_t^2 \\ -\theta_z z_t \\ \end{array} \right]dt + \left[ \begin{array}{c} \sigma_y & 0 \\ 0 & \sigma_z \\ \end{array}\right] dW_t$
where $W_t$ is a Wiener process.
If it simplifies things, what I am ultimately hoping to deduce is the conditional distribution, $y_{t+\Delta t} | z_{t + \Delta t}, y_t$.
I am assuming $dW_{t}=(dW_{t}^{1},dW^{2}_{t})$. This is basically solving a linear process. First we solve the Ornstein–Uhlenbeck process $z_t$ equation to get
$$\large z_{t}= \frac{\sigma_{z}}{\sqrt{2\theta_{z}}} e^{-\theta_{z} t} W^{2}(e^{2 \theta_{z} t})$$
(see https://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_process#Properties_of_sample_paths).
Now you are just left with a linear SDE
$$dy_{t}=(ay_{t}+b(t))dt+\sigma_{y}dW^{1}(t),$$
for $b(t):=\theta_y \alpha + \zeta_0 z_t + \zeta_1 z_t^2 $ and $a:=-\theta_{y}\alpha$, which can be solved as here Solution to General Linear SDE
\begin{align*} Y_t = & Y_0 e^{ at}+ e^{ at}\left( \int_0^t e^{ -as}b(s) \mathrm{d}s + \int_0^t e^{ -as}\sigma_{y} \mathrm{d}B_s\right). \end{align*}