Suppose that under the risk-neutral measure $\mathbf{Q}$ we have an HJM framework dynamics for the instantaneous forward rate $$df_{t,T} = \left(\ldots\right) dt + {}^t \sigma_f (t,T) d W^{Q}_t$$ where $\sigma_f \in\mathbf{R}^d$ is deterministic and where $W^{\mathbf{Q}} \in\mathbf{R}^d$ is a Brownian Motion.
If we suppose that $\sigma_f$ is separable, that is, that you can write $$\sigma_f (t,T) = g(t) h(T)$$ for $g$ a $d\times d$ deterministic matrix function and $h$ a $d$-dimensional deterministic vector function, you can develop a whole theory. Namely, if you set $$H(t) = diag (h(t))$$ and $$\chi (t) = - H'(t) H(t)^1$$ (assuming differentiability and invertibility) you can show that $r_t = f_{0,t} + x_{1,t} + \ldots x_{d,t}$ where the $x_i$ are the coefficient of $d$-dimensional random vector $x_t$ satysfying the following SDE : $$ dx_t = (y(t)\mathbf{1} - \chi(t) x_t) dt + {}^t \sigma_x (t) dW^{\mathbf{Q}}_t$$ with $x_0 = 0$ where $\mathbf{1}$ is the $d$-dimensional vector with all coefficients equal to $1$, $\sigma_x (t) \equiv g(t) h(t)$ and where $y(t)$ is a $d\times d$ deterministic matrix equal to $$H(t) \left( \int_0^t {}^t g(s) g(s) \right) H(t).$$
One can also show that $y(t)$ satisfies the following ODE : $$y'(t) = H(t) {}^t g(t) g(t) H(t) - \chi (t) y(t) - y(t) \chi(t)$$ with $y(0) = 0$.
Now if we want to simulate $x$ at a time discretization $0 = t_0 < t_1 < \ldots < t_d$ one can tactically take advantage of the fact that $x_{t_{i+1}}$ is, conditionally to $x_{t_i}$, a gaussian vector with computable mean and variance, namely : $$(E)\;\;\;\;\;\;\;\;\mathbf{E}^{\mathbf{Q}}\left[ x_{t_{i+1}} \left| x_{t_i} \right.\right] = e^{-\int_{t_i}^{t_{i+1}} \chi (u)du} x_{t_i} + \int_{t_i}^{t_{i+1}} e^{-\int_s^{t_{i+1}} \chi (u)du} y(s) \mathbf{1} ds.$$
and write that $$x_{t_{i+1}} = e^{-\int_{t_i}^{t_{i+1}} \chi (u)du} x_{t_i} + \int_{t_i}^{t_{i+1}} e^{-\int_s^{t_{i+1}} \chi (u)du} y(s) \mathbf{1} ds + \sqrt{\mathbf{Var}^{\mathbf{Q}}\left[ x_{t_{i+1}} \left| x_{t_i} \right.\right]} Z_i$$ where $\sqrt{\mathbf{Var}^{\mathbf{Q}}\left[ x_{t_{i+1}} \left| x_{t_i} \right.\right]}$ is a square-root of the variance matrix (Cholesky for instance) and $Z_1,\ldots,Z_d$ a sequence of independent and identically distributed standard $d$-dimensional gaussian vectors.
Fine. So we need to treat numerically the integral $$\int_{t_i}^{t_{i+1}} e^{-\int_s^{t_{i+1}} \chi (u)du} y(s) \mathbf{1} ds$$ from equation (E).
How do we do that under the hypothesis that all deterministic functions are piecewise constant on the given discretization ? Do we explicitely compute $y(s)$ on $[t_i, t_{i+1}]$ from $y$'s explicit formula or do we use the ODE satisfied by $y$ and if so, how ? Or do we simply say that $$\int_{t_i}^{t_{i+1}} e^{-\int_s^{t_{i+1}} \chi (u)du} y(s) \mathbf{1} ds \simeq (t_{i+1} - t_i) e^{-\int_{t_i}^{t_{i+1}} \chi (u)du} y(t_i) \mathbf{1}$$ and calculate recursively the $y(t_i)$'s from a discretized version of the ODE satisfied by $y$ ?
