I am trying to compute a forward equation of an Ornstein-Uhlenbeck process defined in the line below to use in MATLAB to make pretty pictures
$$dX_{t} = \theta(X_{t}-\mu)dt + \sigma dW_{t}$$
by solving the differential equation using variation of parameters on
$$dX_{t+\Delta t} = \theta(X_{t+\Delta t}-\mu) d(t+h) + \sigma dW_{t+h}.$$
After some work, I end up with the following:
$$X_{t+\Delta t} = e^{\theta \Delta t}X_{t} + \mu(1-e^{\theta\Delta t}) + \sigma e^{\theta(t+\Delta t)}\int_{t}^{t+\Delta t}e^{-\theta s}dW_{s}.$$
From what I've seen online, everything seems fine except the last integral is supposed to look something like
$$\sigma e^{\theta(t+\Delta t)}\int_{t}^{t+\Delta t}e^{-\theta s}dW_{s} = \sigma \sqrt{-\frac{1-e^{-2\theta \Delta t}}{2\theta}}Z_{t},\quad Z_{t}\sim N(0,1)$$
but I don't know how to integrate this with respect to $W_{s}$. I do not know a lot about stochastic calculus and this isn't an assignment so I'm really struggling trying to figure out what tool I need to do to finish this calculation.
One observation I made is that the result I am supposed to get looks like the square root of the result of applying the Ito isometry to the $\int \cdot \,dW_{s}$ term, but I don't understand where the $Z_{t}$ would come from or why the Ito isometry is even relevant to my task if I am not trying to compute the expectation or variance.
My other observation is that if this were a non-stochastic integral, at this point I would probably want to approximate the integral with Gaussian Quadrature or something, but I'm not sure if that is the right way to think of it either.