Consider a matrix differential equation for the vector $\mathbf{y}$ of the form \begin{equation} \dot{\mathbf{y}}(t) = A \, \mathbf{y}(t) + \mathbf{b}(t) \, , \end{equation} where $A$ is a constant square matrix and $\mathbf{b}(t)$ a given time-dependent vector of the same dimension of $\mathbf{y}$. Now, assume that $\mathbf{b}(t)$ is a "white noise" term such that $$ \langle b_i(t) \rangle=0 \qquad \langle b_i(t) b_j(t') \rangle = 2 D_{ij} \, \delta(t-t') $$
where $D$ is a square, positive definite, "diffusion matrix". This process is a multivariate Ornstein–Uhlenbeck process. Is there some good reference for this kind of problem?
We can formally solve the problem by means of the standard matrix exponentiation, namely \begin{equation} {\mathbf{y}(t) } = e^{A t} \mathbf{y}(0) + e^{A t} \int_0^t e^{-A z} \mathbf{b}(z) \, dz \end{equation} Clearly, if $A=0$, we just have the usual Wiener (aka "Brownian") process: \begin{equation} {\mathbf{y}(t) } = \mathbf{y}(0) + \int_0^t \mathbf{b}(z) \, dz \end{equation} and the integral gives us a random variable at time $t$ that is Gaussian with a variance that grows linearly in time $t$. However, I do not know how to deal with the generic integral when $A\neq 0$. References/hints/ideas?
NOTE: I wrote the process in the "Langevin" form. It is immediate to interpret it in the Ito's way by saying that $\mathbf{b}dt = d\mathbf{W}$, where $d\mathbf{W}$ is the increment of the associated multivariate Wiener process (i.e. the "white noise" is the formal time derivative of the Wiener process).
Thanks to the comments I found many references: I added an answer to "close the case" and to list them for the interested ones. This tutorial on stochastic calculus ("Tutorial on Stochastic Differential Equations" by J. R. Movellan) has been proved to be very useful, clear and to the point.
After I was directed in the right direction by the useful comments of @Tobsn, I found nice answers here and here, see also this interesting related question.
Solution in the 1D case: this case is simple and can be immediately understood in terms of Ito's isometry. Direct application of this general result gives exactly the solution proposed by @Tobsn in the comments.
Solution for the multivariate case: thanks to the comment of @kevinkayaks, I found the full solution to all my doubts in section "4.4.6 Multivariate Ornstein-Uhlenbeck Process" of "Handbook of Stochastic Methods" by Gardiner. I checked many references and this stands out for clarity.
Alternatively, "Alternative way to derive the distribution of the multivariate Ornstein–Uhlenbeck process" is about the PDF of the Multivariate Ornstein-Uhlenbeck process. An even better, and very clear, reference is "Fast Bayesian inference of the multivariate Ornstein-Uhlenbeck process".