Suppose we have a system of coupled differential equations in the following form- $$\begin{pmatrix} \dot a\\ \dot b \\\dot c \\ \dot d \end{pmatrix}=\begin{pmatrix} M_{11} & M_{12} & M_{13}& M_{14}\\M_{21} &M_{22} &M_{23} & M_{24}\\ M_{31} &M_{32}& M_{33}&M_{34} \\M_{41} & M_{42} & M_{43} & M_{44} \end{pmatrix} \begin{pmatrix} a \\ b \\ c \\ d\end{pmatrix}+\begin{pmatrix} p \\ q \\r \\s \end{pmatrix}$$
where the matrix $[M_{ij}]_{4 \times 4}$ is time dependent, as well as the matrix of $p,q,r$ and $s$. Suppose we want to solve this set differential equations in the frequency domain instead of the time domain. How can we take its fourier transform? I tried loking up DFT but couldn't get any satisfactory answers on this.
We can write the matrix differential equation as $\dot{\mathbf{u}} = \mathbf{M} \mathbf{u} + \mathbf{c}$ where $$ \mathbf{u} = \begin{pmatrix}a\\b\\c\\d\end{pmatrix},\quad \mathbf{M} = \begin{pmatrix}M_{11} & M_{12} & M_{13}& M_{14}\\M_{21} &M_{22} &M_{23} & M_{24}\\ M_{31} &M_{32}& M_{33}&M_{34} \\M_{41} & M_{42} & M_{43} & M_{44}\end{pmatrix},\quad \mathbf{c} = \begin{pmatrix}r\\q\\r\\s\end{pmatrix} . $$
Taking the Fourier transform of both sides results in $i\omega\,\hat{\mathbf{u}} = \hat{\mathbf{M}}*\hat{\mathbf{u}} + \hat{\mathbf{c}}$ where $\hat{}$ denotes Fourier transform and $*$ denotes convolution. Here the Fourier transforms of the matrices are done component-wise, i.e. $$ \hat{\mathbf{u}} = \begin{pmatrix}\hat a\\\hat b\\\hat c\\\hat d\end{pmatrix},\quad \hat{\mathbf{M}} = \begin{pmatrix}\hat M_{11} & \hat M_{12} & \hat M_{13}& \hat M_{14}\\\hat M_{21} &\hat M_{22} &\hat M_{23} & \hat M_{24}\\ \hat M_{31} &\hat M_{32}& \hat M_{33}&\hat M_{34} \\\hat M_{41} & \hat M_{42} & \hat M_{43} & \hat M_{44}\end{pmatrix},\quad \hat{\mathbf{c}} = \begin{pmatrix}\hat r\\\hat q\\\hat r\\\hat s\end{pmatrix} . $$
Thus you get a non-differential equation, but with a convolution which makes things a bit difficult.