Let $A:=\text{diag}(a_1,\dots,a_N)$, the diagonal matrix with entries $a_1,\dots,a_N$. Let $U(t):=B\circ W(t)$, where $B:=(b_{ij})$ (the matrix of the $b_{ij}$'s), $W(t):=(W_{ij}(t))$ (the standard Wiener process in $\mathbb{R}^{N\times N}$), and $\circ$ is the Hadamard matrix product, so that $U_{ij}(t)=b_{ij}W_{ij}(t)$ for all $i,j,t$.
Then it is rather straightforward to show that the quadratic variation of $U$ is
\begin{equation}
[U]_t=\lim\sum_{r=0}^{n-1}(U(t_{r+1})-U(t_r))^2=t\tilde B^2
\tag{1}
\end{equation}
for $t>0$,
where $0=t_0<\dots<t_n=t$, the limit is in $L^2$ (and hence in probability) as $\max_r(t_{r+1}-t_r)\to0$, and
$$\tilde B^2:=\text{diag}(b_{11}^2,\dots,b_{NN}^2).$$
In general, there is apparently no explicit solution of the system of equation $dX_i=X_i a_i \, dt+\sum_{j=1}^N X_j b_{ij} \,dW_{ij}$; cf. the last line on page 2 of Duan--Yan.
Consider the special case when $A=aI_N$ and $\tilde B^2=b^2I_N$ for some real $a$ and $b$, and let
\begin{equation}
X(t):=e^{(a-b^2/2)t}e^{B\circ W(t)}X_0 \tag{2}
\end{equation}
with $X_0$ independent of $W(\cdot)$. Then it follows from (1) (say, along the lines of the derivation of the Ito formula for real-valued functions of a real-valued Ito process) that (2) is a solution to the system of equations $dX_i=X_i a_i \, dt+\sum_{j=1}^N X_j b_{ij} \,dW_{ij}$.
The case of the displayed system of two stochastic differential equations in question is covered by this setting, with $N=2$, $a=0$, and $b=1$.
