I need to write down a code to simulate an SDE of the following type: $$ dX_t = - ( \eta_t X_t + \chi_t ) dt + \sigma dW_t, \ X_0=v_0 \quad t \in [0,T] $$ where $\eta_t$ is a deterministic function, let us say $\eta_t = t, t \in [0,T]$,$X_t \in \mathbb{R}^d$, $\sigma, v_0$ are two scalars, and $\chi_t$ is of the form $$ \chi_t^i = \int_t^T e^{\int_s^t \eta_u du}ds, \quad t \in [0,T], i =1 , \dots , d. $$ I also know that is possible to give a closed formula for the explicit solution which is $$ X_t = e^{- \int_0^t \eta_u du} v_0 - \int_0^t e^{-\int_s^t \eta_u du}\chi_s ds + \sigma \int_0^t e^{- \int_s^t\eta_ud_u} dW_s. $$ Just to be more specific with the physical meaning of $X_t$ it represents a velocity, hence I am also interested in numerically simulate $x_0 +\int_0^t X_s ds$ which is the position at time $t$.
Given all these not so clear premises, I don't know anything about simulating SDE and I was wondering if someone can help me in finding some reference to simulate such a problem with matlab.
To deep dive into numerical methods for SDEs, a good start will be the Euler-Mayurama method to simulate the solution of the SDE. The headlines of this method are the following: Denoting $X_t = \lbrace{X_t^j\rbrace}_{j\in\lbrace{1,...,d}\rbrace}$, we have \begin{align} X^j_{t_{i+1}} = X^j_{t_{i}} - (\eta_{t_{i}}X^j_{t_{i}}+\chi_{t_{i}})\Delta_{t_i} + \sigma(W^j_{t_{i+1}}-W^j_{t_{i}})\\ = X^j_{t_{i}} - (\eta_{t_{i}}X^j_{t_{i}}+\chi_{t_{i}})\Delta_{t_i} + \sigma\sqrt{t_{i+1}-t_i}\epsilon^j \end{align} Where $X_{t_0} = X_0$, $\Delta_{t_i} = t_{i+1}-t_{i}$ for $i = 0, ..., N$ and $\epsilon$ is a Gaussian vector with mean $0$ and co-variance $I_d$. It is an iterative method as the solution of the SDE is updated at each time step.
Another way will be to use the explicit solution as given in your post. The integral $\int_0^t\exp(-\int_s^t\eta_udu)\chi_sds$ can be computed using well-known numerical methods and the stochastic integral, i.e.$\left(\sigma\int_0^t\exp(-\int_s^t\eta_udu)dW_s\right)^j$, can be computed using the definition of the stochastic integral.
Once you have simulated $X_t$ using the methods described above, computing $x_0 +\int_0^t X_s ds$ could be done as follows: \begin{equation} x_0 +\int_0^t X_s ds = x_0 + \sum_{i=0}^N X_{t_i}\Delta_{t_i} \end{equation}
Note that the Euler-Mayurama method is well described in this document from Giles Pages (Chapter 7). You can also find many numerical methods applied in the probability field in the same doc.