I am trying simulate a to solve a 2-dimensional stochastic process and $Y_t^1$ is a mean-reverting square root process which I simulated on a time grid using its known conditional distribution. I want to compute the stochastic process $Y_t^2$ given below . \begin{eqnarray*} Y_t^2:=e^{\tilde{a}_{22}t}\tilde{s}_2+\int_0^te^{\tilde{a}_{22}(t-s)}\left[\left(\tilde{a}_{21}Y_s^1+ \tilde{b}_2\right)ds+\sqrt{ \alpha+\vert \beta\vert Y_s^1}d\tilde{W}^2_s\right] \end{eqnarray*}
That means that I want to compute the values of $Y_t^2$ along 100 paths at the time discretization grid (we know that the time discretization grid is $[0 ,0.01, 0.02 ,\dots,0.99]$. I want to compute the a $100 \times 100$ matrix where the columns vector represent the paths. $\tilde{a}_{21}, \tilde{a}_{22},\alpha, \vert \beta\vert,\tilde{b}_2,\tilde{s}_2$ are known constants and $Y_t^1$ is known (in other words $100 \times 100$ matrix with each column representing the paths of $Y_t^1$) is known.
Attempt: I tried to approximate the integral by using the following approximation with
$t_1=0,t_2=0.01,t_3=0.02,\dots,t_{100}=0.99$
$Y_{t_k}^2(\omega)\approx e^{\tilde{a}_{22}t_k}\tilde{s}_2+\tilde{a}_{21} Y_0^1(\omega)\int_0^{t_2} e^{\tilde{a}_{22}(t_k-s)}ds+\tilde{b}_2\int_0^{t_2}e^{\tilde{a}_{22}(t_k-s)}ds+\sqrt{\alpha+\vert \beta \vert Y_0^1(\omega)}\int_0^{t_2}e^{\tilde{a}_{22}(t_k-s)}d\tilde{W}^2_s+\dots +\tilde{a}_{21} Y_{t_{k-1}}^1(\omega)\int_{t_{k-1}}^{t_k} e^{\tilde{a}_{22}(t_k-s)}ds+\tilde{b}_2\int_{t_{k-1}}^{t_k} e^{\tilde{a}_{22}(t_k-s)}ds+\sqrt{\alpha+\vert \beta \vert Y_{t_{k-1}}^1(\omega)}\int_{t_{k-1}}^{t_k}e^{\tilde{a}_{22}(t_k-s)}d\tilde{W}^2_s$.
I compute the stochastic integral by assuming the integrand takes value at the left-endpoint of the interval which is $t_1=0$ in the case below and integrand comes out of the integral and we are just left with stochastic integral the Brownian motion which is simulated by simulating normally distributed random variables of mean zero and variance equal to the time step(in this case $0.01$)
$ \sqrt{\alpha+\vert \beta \vert Y_0^1(\omega)}\int_0^{t_2}e^{\tilde{a}_{22}(t_k-s)}d\tilde{W}^2_s=\sqrt{\alpha+\vert \beta \vert Y_0^1(\omega)}e^{\tilde{a}_{22}t_k}\int_0^{t_2}d\tilde{W}^2_s\\=\sqrt{\alpha+\vert \beta \vert Y_0^1(\omega)}e^{\tilde{a}_{22}t_k} *N(0,0.01) \text{sample} $
I compute the lebesgue integral as follows
$\tilde{a}_{21} Y_{t_{k-1}}^1(\omega)\int_{t_{k-1}}^{t_k} e^{\tilde{a}_{22}(t_k-s)}ds= Y_{t_{k-1}}^1(\omega)\frac{\tilde{a}_{21}}{\tilde{a}_{22}}\left(e^{\tilde{a}_{22}(t_k-t_{k-1})}-e^{\tilde{a}_{22}(t_k-t_k)}\right)$