A question about Finite Element Method. Can we solve a matrix differential equation $\textbf{X}'(t)=\textbf{A}\textbf{X}(t)+\textbf{B}(t)$ by Finite Element Method?
I will be very happy if you give some special example or suggest some books about the method for solving matrix differential eauations.
In here, 
\begin{equation}
\begin{aligned}
&\textbf{X}(t)=\left( \begin{array}{cccc}
x_{1}(t) \\
x_{2}(t) \\
x_{3}(t)\\
x_{4}(t) \end{array} \right), \textbf{A}
=
\left( \begin{array}{cccc}
a_{11} & a_{12} &a_{13} & a_{14}\\
a_{21} & a_{22} & a_{23} & a_{24} \\
a_{31} & a_{32}& a_{33} & a_{34} \\
a_{41} & a_{42} & a_{43} & a_{44}
\end{array} \right),\\
&\textbf{B}(t)=\left( \begin{array}{cccc}
b_{1}(t)\\
b_{2}(t)\\
b_{3}(t)\\
b_{4}(t) \end{array} \right).
\end{aligned}
\end{equation}
initial conditions:
$$x_1(0)=m_1$$
$$x_2(0)=m_2$$
$$x_3(0)=m_3$$
$$x_4(0)=m_4.$$
 A: The differential system ${\bf X}' = {\bf A} {\bf X} + {\bf B}$ where $t \mapsto {\bf X}(t)$ is a vector-valued function of size $n=4$, ${\bf A}$ is a $n$-by-$n$ matrix and $t \mapsto {\bf B}(t)$ is a vector-valued function of size $n$, can be solved by hand in many particular cases. Indeed, this system is an order-1 linear system of ordinary differential equations (ODE). Solutions write as
$$
{\bf X}(t) = \exp(t{\bf A})\, {\bf X}(0) + \int_{0}^t \exp((t-\tau) {\bf A}) \, {\bf B}(\tau)\, \mathrm{d}\tau \, ,
$$
where $\exp(s{\bf A})$ is a matrix exponential. 
Therefore, if the matrix exponential and the integral can be computed, then the solution is given above.
(1) FEM. To implement FEM, the first step is to write a weak formulation of the ODE. We write the scalar product of the ODE with the test vector field ${\bf Y}(t)$, and we integrate over $t\in \Theta$. Hence, using integration by parts,
$$
\int_\Theta ({\bf Y}' +  {\bf A}^\top{\bf Y})^\top\, {\bf X} = -\int_\Theta {\bf Y}^\top {\bf B},
$$
for all ${\bf Y}$ which trace equals zero on the boundaries of $\Theta$. One notes that the above weak formulation writes as $a({\bf X},{\bf Y}) = f({\bf Y})$, where $a(\cdot,\cdot)$ is a bilinear form over $H^1(\Theta)$ and $f(\cdot )$ is a linear form. It remains to discretize the above equations by considering that ${\bf X}$, ${\bf Y}$ belong to finite-dimension spaces, e.g. by using appropriate bases of Lagrange polynomials of the variable $t$.
(2) FD. The system ${\bf X}' = {\bf A} {\bf X} + {\bf B}$ can be integrated numerically by finite-difference methods. The simplest method is the first-order accurate Euler method. Here, the forward Euler method writes
$$
\left\lbrace
\begin{aligned}
&{\bf X}_{n+1} = {\bf X}_{n} + \Delta t \left({\bf A}{\bf X}_{n} + {\bf B}(t_n)\right) , \\
&{\bf X}_{0} = {\bf X}(0) \, ,
\end{aligned}
\right.
$$
where ${\bf X}_{n} \simeq {\bf X}(t_{n})$, and $\Delta t = t_{n+1} - t_n$ denotes the time step. The method is stable provided that
$
\Delta t \leq {2}/{\varrho({\bf A})}
$,
where $\varrho({\bf A})$ is the spectral radius of the matrix $\bf A$. The backward Euler method
$$
\left\lbrace
\begin{aligned}
&{\bf X}_{n+1} = \left({\bf I} - \Delta t {\bf A}\right)^{-1}\, \left({\bf X}_{n} + \Delta t\, {\bf B}(t_{n+1})\right) , \\
&{\bf X}_{0} = {\bf X}(0) \, ,
\end{aligned}
\right.
$$
where ${\bf I}$ is the identity matrix, is unconditionally stable. Thus, it bypasses stability issues, but requires the computation of a matrix inverse.
Many more accurate numerical methods exist, such as Runge-Kutta methods. If the system is stiff, then the time step $\Delta t$ must be chosen very small to get good accuracy, which is costly from a computational point of view. Hence, dedicated adaptive methods may be preferred. For instance, the reader is referred to the book Numerical Recipes in C: The Art of Scientific Computing for details.
A: This is perfectly possible. First we need a weak formulation: take smooth $v$ then 
$$
\int_0^T v^T(x' - Ax) dt = \int_0^T v^T Bdt
$$
is a nice weak formulation. It remains to choose the spaces $X_h$ and $V_h$ to take $x$ and $v$ from.
If $X_h=V_h$=piecewise polynomials, then the resulting system completely couples over all time steps, which is not nice.
A sensible choice is to take one of $X_h$ and $V_h$ piecewise polynomial and continuous, the other space piecewise polynomial but discontinuous. You end up with a Petrov-Galerkin or discontinuous Galerkin scheme. Due to the discontinuity in one of the variables, the time steps are decoupled. 
Choosing piecewise linear and continuous functions for $x$, but discontinuous and piecewise constant functions for $v$ results in a variant of the Crank-Nicolson scheme.
