\documentclass[reqno]{amsart}
\usepackage{hyperref}
\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2008(2008), No. 95, pp. 1--12.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu (login: ftp)}
\thanks{\copyright 2008 Texas State University - San Marcos.}
\vspace{9mm}}
\begin{document}
\title[\hfilneg EJDE-2008/95\hfil Linear matrix differential equations]
{Linear matrix differential equations of higher-order and
applications}
\author[R. Ben Taher, M. Rachidi \hfil EJDE-2008/95\hfilneg]
{Rajae Ben Taher, Mustapha Rachidi} % in alphabetical order
\address{Rajae Ben Taher \newline
D\'epartement de Math\'ematiques et Informatique \\
Facult\'e des Sciences, \newline
Universit\'e Moulay Ismail \\
B.P. 4010, Beni M'hamed, M\'ekn\'es, Morocco}
\email{bentaher@fs-umi.ac.ma}
\address{Mustapha Rachidi \newline
Mathematics Section LEGT - F. Arago, Acad\'emie de Reims \\
1, Rue F. Arago, 51100 Reims, France}
\email{mu.rachidi@hotmail.fr}
\thanks{Submitted February 25, 2007 Published July 5, 2008.}
\subjclass[2000]{15A99, 40A05, 40A25, 15A18}
\keywords{Algebra of square matrices;
linear recursive relations; matrix powers; \hfill\break\indent
matrix exponential; linear matrix differential equations}
\begin{abstract}
In this article, we study linear differential equations
of higher-order whose coefficients are square matrices.
The combinatorial method for computing the matrix powers
and exponential is adopted. New formulas representing
auxiliary results are obtained. This allows us to prove
properties of a large class of linear matrix differential
equations of higher-order, in particular results of Apostol and
Kolodner are recovered. Also illustrative examples and
applications are presented.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{example}[theorem]{Example}
\newtheorem{remark}[theorem]{Remark}
\section{Introduction}\label{section1}
Linear matrix differential equations are important in many fields of
mathematics and their applications. Among the most simple and
fundamental are the first order linear systems with constant
coefficients,
\begin{equation}\label{Equadiff-12}
X'(t)=AX(t), \quad \text{such that } X(0)=1_d,
\end{equation}
where $1_d$ is the identity matrix of $\mathcal{{M}}(d; \mathbb{C})$,
$A\in \mathcal{M}(d; \mathbb{C})$ and
$X\in {\mathcal{C}}^{\infty}(\mathbb{R}; \mathcal{{M}}(d; \mathbb{C}))$.
For the sake of simplicity we set in
the sequel ${ \mathcal{A}}_d=\mathcal{{M}}(d; \mathbb{C})$. The system
(\ref{Equadiff-12}) has been extensively studied;
its solutions depend closely on the computation of
$e^{tA}$ ($t\in \mathbb{R}$). To perform this computation, many
theoretical and numerical methods have been developed (see
\cite{br1, br2, br3, cy, ga, leo, mv, mr2, rp, vs} for example).
The combinatorial method is among those considered recently for
obtaining some practical expressions of $A^n$ ($n\geq d$) and
$e^{tA}$ ($t\in \mathbb{R}$) (see \cite{br1, br2, mr2}). Techniques of
this method are based on properties of some linear recursive
sequences (see \cite{le, mr1} for example), known in the
literature as $r$-generalized Fibonacci sequences (to abbreviate
we write $r$-GFS).
Let $A_0,\dots, A_{s-1}$ be in ${ \mathcal{A}}_d$ and consider
the linear matrix differential equation of higher-order
\begin{equation}\label{Equadiff-1}
X^{(s)}(t)=A_0X^{(s-1)}(t)+\dots+A_{s-1}X(t),
\end{equation}
subject to the initial conditions $X(0)$, $X'(0)$,
$\dots$, $X^{(s-1)}(0)$. To study \eqref{Equadiff-1} it is customary
to write it under its
equivalent form as the linear first order system
(\ref{Equadiff-12}). More precisely, \eqref{Equadiff-1}
takes the form $Y'(t)=BY(t)$, where $B\in {\mathcal{A}}_{ds}$ and
$Y \in {\mathcal{C}}^{\infty}(\mathbb{R};{\mathcal{A}}_{ds})$ (see Section
4 for more details). For $s=2$ properties of the linear
matrix differential equation \eqref{Equadiff-1},
\begin{equation}\label{Equadiff-3}
X''(t)=A_0X'(t)+A_1X(t),
\end{equation}
have been studied (directly) by Apostol and Kolodner
(see \cite{ap, ko}) without
appealing to its equivalent form as equation (\ref{Equadiff-12}).
The main purpose of this paper is to study solutions of a large
class of linear matrix differential equations \eqref{Equadiff-1},
whose coefficients are in the algebra $\mathcal{A}_d$. First
we consider \eqref{Equadiff-1} when
$A_0=\dots=A_{r-2}=\Theta_d$ (the zero matrix of
$\mathcal{A}_d$) and $A_{r-1}=A\neq \Theta_d$. Solutions to
these matrix differential equations are expressed in terms
of the coefficients of the polynomial $P(z)=z^r-a_0z^{r-1}-\dots
-a_{r-1}$ satisfying $P(A)=\Theta_d$ and matrices $A^0=1_d$, $A$,
$\dots$, $A^{r-1}$. Moreover, these solutions are also described
with the aid of the zeros of the polynomial $P(z)$. Furthermore,
the case $r=2$ is improved and the fundamental results of
Apostol-Kolodner are recovered and their extension is established.
Secondly, in light of a survey of the general equation
\eqref{Equadiff-1}, we manage to supply their solutions under a
combinatorial form.
In Section 2 we recall auxiliary
results on the powers and exponential of an element
$A\in\mathcal{A}_d$. In Section 3, we study the class
of linear matrix differential equations of higher-order
\eqref{Equadiff-1} associated to $A_0=\dots=A_{s-2}=\Theta_d$ and
$A_{s-1}=A$ (we call this class the {\it Higher order linear
matrix differential of Apostol type}) and recover the results of
Apostol-Kolodner. In Section 4 we consider the combinatorial aspect
of the matrix powers and exponential in order to
explore the general setting of linear matrix differential
equations \eqref{Equadiff-1}. Finally,
concluding remarks are stated and a future problem is formulated
in Section 5.
\section{Auxiliary results on the powers and exponential of a matrix}
\label{section2}
We review here some basic results on the matrix powers and
exponential needed in the next sections. That is, we recall some
results of \cite{br1, br2, br3} and set forth a new result in
Proposition \ref{prop-24}. To begin, let $A_0,\dots, A_{r-1}$ and
$S_0,\dots, S_{r-1}$ be two finite sequences of $\mathcal{A}_d$
with $A_{r-1}\neq \Theta_d$. Consider the recursive sequence
$\{Y_n\}_{n\geq 0}$ such that $Y_n=S_n$ for $0\leq n\leq r-1$ and
\begin{equation}\label{recursive-1}
Y_{n+1}=A_0Y_n+\dots+A_{r-1}Y_{n-r+1}, \quad \text{for every } n\geq r-1.
\end{equation}
When $A_jA_k=A_kA_j$, for every $0\leq j,k\leq r-1$ and
following the same straightforward computation as in \cite{br1, mr2}
we obtain,
\begin{equation}\label{combinat-1}
Y_n=\rho(n,r)W_0+\dots+\rho(n-r+1,r)W_{r-1}, \quad\text{for every }
n\geq r,
\end{equation}
where $W_p=A_{r-1}S_p+\dots+A_pS_{r-1}$ ($p=0,1,\dots, r-1$) and
\begin{equation}\label{rho-1}
\rho(n,r)=\sum_{k_0+2k_1+\dots +rk_{r-1}=n-r} \frac{(k_0+\dots
+k_{r-1})!}{k_0!k_1!\dots k_{r-1}!}A_0^{k_0}A_1^{k_1}\dots
A_{r-1}^{k_{r-1}},
\end{equation}
for every $n\geq r$, with $\rho(r,r)=1$ and $\rho(n,r)=0$
for $n\leq r-1$ (see \cite{br1, br2, leo, mr1, mr2}). The
computation of the powers $A^n$ ($n\geq r$) follows by a direct
application of \eqref{combinat-1}-\eqref{rho-1}. Indeed, the
polynomial equation $P(A)=\Theta_d$ shows that $A^{n+1}=
a_0A^{n}+\dots+a_{r-1}A^{n-r+1}$, for every $n\geq r-1$, and by the
way the sequence $\{A^n\}_{n\geq 0}$ is an $r$-GFS in
$\mathcal{A}_d$, whose coefficients and initial values are
$A_0=a_0,\dots , A_{r-1}=a_{r-1}$, $S_0=A^0=1_d, \dots,
S_{r-1}=A^{r-1}$ (respectively). In light of \cite{br1, mr2} we have
the following characterization of the powers of $A$,
\begin{equation}\label{powers-1}
A^n=\sum_{p=0}^{r-1}\Big(\sum_{j=0}^pa_{r-p+j-1}\rho(n-j,r)\Big)A^p,\;
\text{for every } n\geq r,
\end{equation}
where
\begin{equation}\label{rho-2}
\rho(n,r)=\sum_{k_0+2k_1+\dots+rk_{r-1}=n-r}\frac{(k_0+\dots+k_{r-1})!}{k_0!k_1!\dots
k_{r-1}!}a_0^{k_0}\dots a_{r-1}^{k_{r-1}}, \quad \text{for }
n\geq r,
\end{equation}
with $\rho(r,r)=1$ and $\rho(n,r)=0$ for $n\leq r-1$
(see \cite{br1, mr2}). The matrix exponential
$e^{tA}$ ($A\in\mathcal{A}_d$) is
defined as usual by the series
$e^{tA}=\sum_{n=0}^{+\infty}(t^n/n!)A^n$. It turns out following
\cite{br1} that direct computation using
(\ref{powers-1})-\eqref{rho-2} allows us to derive that the unique
solution of the matrix differential equation (\ref{Equadiff-12}),
satisfying the initial condition $X(0)=1_d$, is
$X(t)=e^{tA}=\sum_{k=0}^{r-1}\Omega_k(t)A^k$ such that
\begin{equation}\label{omega-1}
\Omega_k(t)=\frac{t^k}{k!}+\sum_{j=0}^ka_{r-k+j-1}\omega_j(t),\quad
(0\leq k\leq r-1),
\end{equation}
with
\[
\omega_j(t)=\sum_{n=0}^{+\infty}\rho(n-j,r)\frac{t^n}{n!}
\quad (0\leq j\leq r-1)
\]
and $\rho(n,r)$ given by \eqref{rho-2}.
As shown here the preceding expressions of $A^n$ ($n\geq r$) and
$e^{tA}$ are formulated in terms of the coefficients of the
polynomial $P(z)=z^r-a_0z^{r-1}-\dots-a_{r-1}$ and the first $r$
powers $A^0=1_d$, $A$, $\dots$, $A^{r-1}$. The goal now is to
express $A^n$ ($n\geq r$) and $e^{tA}$ using the roots of the
polynomial $P(z)=z^r-a_0z^{r-1}-\dots -a_{r-1}$ satisfying
$P(A)=\Theta_d$. To this aim, let $\lambda_1, \dots , \lambda_s$
be the distinct roots of $P(z)=z^r-a_0z^{r-1}-\dots-a_{r-1}$ of
multiplicities $m_1,\dots,m_{s}$ (respectively). For every $j$
($1\leq j\leq s$) we consider the sequence $\{b_{i,j}\}_{0\leq
i\leq m_j-1}$ such that $b_{i,j}=0$ for
$i>m_1+\dots+m_{j-1}+m_{j+1}+\dots+m_s$ and otherwise
\begin{equation}\label{bij}
b_{i,j}=\sum_{h_1+\dots+h_{j-1}+h_{j+1}+\dots+h_s=i,h_d\leq m_d}
\prod_{d=1,d\neq j}^s(^{h_d}_{m_d})(\lambda_j-\lambda_d)^{m_d-h_d}.
\end{equation}
For every $j$ ($1\leq j\leq s$), let $\{\alpha_{i,j}\}_{0\leq i\leq
m_j-1}$ be the sequence defined by $\alpha_{0,j}=1$ and
\begin{equation}\label{alpha-ij}
\alpha_{i,j}=\frac{-1}{b_{0,j}}
\left(b_{1,j}\alpha_{i-1,j}+b_{2,j}\alpha_{i-2,j}+\dots+b_{i-1,j}
\alpha_{1,j}+b_{i,j}\alpha_{0,j}\right),
\end{equation}
where the $b_{i,j}$ are given by $(\ref{bij})$. Recall that
$\{b_{i,j}\}_{0\leq i\leq m_j-1}$ and $\{\alpha_{i,j}\}_{0\leq i\leq
m_j-1}$ have been introduced in \cite{br2} for computing the powers
of a matrix $A\in \mathcal{A}_d$. Besides the exploration of
Expression (4.15) in \cite{vs1} and Proposition 4.3 in \cite{br3}
allow us to obtain an explicit formula for the $\alpha_{i,j}$ as
follows,
\begin{equation}\label{alpha-3}
\alpha_{i,j}=(-1)^i\sum_{h_1+h_2
+ \dots +h_{j-1}+ h_{j+1}\dots+h_s=i,h_d \leq
m_d}\prod_{{d=1}_{d\neq j}}^s(^{h_d}_{m_d+h_d-1})\left(\lambda_j-\lambda_d
\right)^{- h_d }.
\end{equation}
With the aid of the basis of the Lagrange-Sylvester interpolation
polynomials as in \cite{ga}, we obtain the following result.
\begin{proposition}\label{prop-23}
Let $A$ be in $\mathcal{A}_d$ such that
$p(A)=\prod_{j=1}^s(A-\lambda_j1_d)^{m_j}=\Theta_d$
($\lambda_i\neq \lambda_j$, for $i\neq j$). Then
\begin{equation}\label{exp-1}
\begin{aligned}
e^{tA}& = \sum_{j=1}^se^{\lambda_jt}\sum_{k=0}^{m_j-1}f_{j,k}(t)
(A-\lambda_j1_d)^k q_j(A), \\
& = \sum_{j=1}^se^{\lambda_jt}\Big[\sum_{k=0}^{m_j-1}\frac{t^k}{k!}
(A-\lambda_j1_d)^k\Big]p_j(A),
\end{aligned}
\end{equation}
\begin{equation}\label{pow-1}
A^n =
\sum_{j=1}^s\sum_{k=0}^{\min(n,m_j-1)}
\sum_{i=0}^{m_j-k-1}\lambda_j^{n-k}\alpha_{i,j}
\binom{n}{k}
(A-\lambda_j1_d)^{k+i} q_j(A)
\end{equation}
such that
\[
p_j(z)=\prod_{d=1, d\neq j}
((z-\lambda_d)^{m_d}/(\lambda_j-\lambda_d)^{m_d})
\sum_{i=0}^{m_j-1}\alpha_{i,j}(z-\lambda_j)^i,
\]
$q_j(z)=p(z)/(z-\lambda_j)^{m_j}$ and
\[
f_{j,k}(t)=\sum_{i=0}^k\alpha_{i,j}(t^{k-i}/(k-i)!)(1/\prod_{d=1,
d\neq j}(\lambda_j-\lambda_d)^{m_d}),
\]
where the $\alpha_{i,j}$ are given by \eqref{alpha-3}.
\end{proposition}
When $p(A)=(A-\lambda_j1)^{r}=\Theta_d$ we can take $\alpha_{01}=1$
and $\alpha_{i1}=0$ for every $i\neq 0$.
Let $\Omega_0(t)$ be the coefficient of $A^0=1_d$ in the well known
polynomial decomposition
$e^{tA}=\sum_{k=0}^{r-1}\Omega_k(t)A^k$. A direct computation, using
(\ref{exp-1}), yields a new formula for $\Omega_0(t)$ in terms of
$\{\lambda_j\}_{\{1\leq j\leq s\}}$ as follows.
\begin{proposition}\label{prop-24}
Under the hypothesis of Proposition \ref{prop-23}, we have
\begin{equation}\label{omega-0}
\Omega_0(t)=\sum_{k=1}^se^{\lambda_kt}Q_k(t)\prod_{j=1,
j\neq k}^s\frac{(-\lambda_j)^{m_j}}{(\lambda_k-\lambda_j)^{m_j}},
\end{equation}
where
\[
{Q_k(t)=\sum_{p=0}^{m_k-1}\Big(\sum_{i=0}^p\alpha_{i,k}t^{p-i}/(p-i)!\Big)
(-\lambda_k)^p},
\]
with the $\alpha_{i,k}$ given by \eqref{alpha-3}. Moreover, we have
${\frac{d\Omega_{k+1}}{dt}(t)=a_{r-k-2}\frac{d\omega_0}{dt}(t)+\Omega_k(t)}$,
with $\Omega_0(t)=1+a_{r-1}\omega_0(t)$.
\end{proposition}
From Proposition \ref{prop-24} we derive that
\[
\omega_0(t)=\frac{1-\Omega_0(t)}{(-\lambda_1)^{m_1}\dots
(-\lambda_s)^{m_s}},
\]
where $\Omega_0(t)$ is given by (\ref{omega-0}).
It seems for us that (\ref{omega-0}), which gives
$\Omega_0(t)$ in terms of the eigenvalues
$\{\lambda_j\}_{\{1\leq j\leq s\}}$ of the matrix $A$, is not known
in the literature under this form.
The preceding results will play a central role in the next sections
devoted to some properties of the linear matrix differential
equations of higher-order \eqref{Equadiff-1}-\eqref{Equadiff-3}.
\section{Higher-order matrix differential equations of Apostol type}
\label{section3}
We are concerned here with the higher-order matrix differential
equations \eqref{Equadiff-1} when $A_0=\dots =A_{r-2}=\Theta_d$
and $A_{r-1}=A\neq \Theta_d$, where $A\in \mathcal{A}_d$. These
linear matrix differential equations are called {\it of Apostol
type}. For reason of clarity we proceed as follows. We start by
reconsidering the case $r=2$ and afterwards we focus on the case
$r\geq 2$, particularly results of Apostol and Kolonder are
recovered.
\subsection{Simple second-order differential equations}
\label{Subsection-3.1}
Let $A$ be in $\mathcal{A}_d$ satisfying the polynomial equation
$P(A)=\Theta_d$, where $P(z)=z^r-a_0z^{r-1}-\dots-a_{r-1}$. It is
well known that the second-order linear matrix differential equation
\begin{equation}\label{Equadiff-S3}
X''(t)=AX(t),
\end{equation}
has a unique solution ${X(t)=C(t)X(0)+S(t)X'(0)}$,
where $X(0)$ and $X'(0)$ are the prescribed initial values and
$C(t)$ and $S(t)$ are the following series
\[
{C(t)=\sum_{k=0}^{+\infty}\frac{t^{2k}}{(2k)!}A^k}, \quad
{S(t)=\sum_{k=0}^{+\infty}\frac{t^{2k+1}}{(2k+1)!}A^k}
\]
(see \cite{ap, ko}). If we substitute for $A^k$ its expression
given in (\ref{powers-1}) we manage to have the following
improvement of the Apostol-Kolodner result.
\begin{proposition}\label{prop-31}
Let $A$ be in $\mathcal{A}_d$ satisfying the polynomial equation
$P(A)=\Theta_d$, where $P(z)=z^r-a_0z^{r-1}-\dots-a_{r-1}$. Then,
the unique solution of the matrix differential equation
\eqref{Equadiff-S3}, with the prescribed initial values $X(0)$ and
$X'(0)$ is $X(t)=\sum_{k=0}^{r-1}(C_k(t)X(0)+ S_k(t)X'(0))A^k$, such
that
$$
C_k(t)=\frac{t^{2k}}{(2k)!}+\sum_{n=r}^{+\infty}\frac{t^{2n}}{(2n)!}
\rho_k(n),\quad
S_k(t)=\frac{t^{2k+1}}{(2k+1)!}+\sum_{n=r}^{+\infty}
\frac{t^{2n+1}}{(2n+1)!}\rho_k(n)
$$
for $0\leq k\leq r-1$, where
${\rho_k(n)=\sum_{j=0}^ka_{r-k+j-1}\rho(n-j,r)}$ with
the $\rho(n,r)$ given by \eqref{rho-2}.
\end{proposition}
Consider the first-order differential equation (\ref{Equadiff-12})
with $A \in \mathcal{A}_d$ satisfying $P(A)=\Theta_d$, where
$P(X)=z^r-a_0z^{r-1}-\dots-a_{r-1}$. Kolodner's method (see
\cite{ko}) shows that $P(D)X(t)=P(A)X(t)=0$, where $D=d/dt$.
Therefore, we have $X(t)=\sum_{j=0}^{r-1}\phi_j(t)A^j$, where the
functions $\phi_j(t)$ ($0\leq j\leq r-1$) verify the scalar
differential equation $P(D)x(t)=0$, whose initial conditions are
$D^k\phi_j(0)=\delta_{j,k}$ (= 1 for $j=k$ and 0 if not), (see
\cite{ko, leo}). For the linear matrix differential equation
\eqref{Equadiff-S3} we have $D^2X(t)=AX(t)$, and the preceding
method shows that $P(D^2)X(t)=0$. If we set $Q(z)=P(z^2)$, it is
easy to show that each function $\phi_j(t)$ ($0\leq j\leq 2r-1$)
is also a solution of the scalar ordinary differential equation
(of order $2r$) $Q(D)x(t)=0$, satisfying the initial conditions
$D^k\phi_j(0)=\delta_{j,k}$ ($0\leq j\leq 2r-1$). Therefore,
Equation (\ref{omega-1}) implies that
\begin{equation}\label{solu-S32}
\phi_k(t)=\frac{t^k}{k!}+\sum_{j=0}^kb_{2r-k+j-1}\omega_j(t),
\end{equation}
where $b_{2i}=0$, $b_{2i+1}=a_i$ and
$\omega_j(t)=\sum_{n=2r}^{+\infty}\rho(n-j,2r)t^n/n!$ ($0\leq j\leq
2r-1$), with $\rho(n,2r)$ given by \eqref{rho-2} are such that the
$a_i$ and $r$ are replaced by the $b_i$ and $2r$ (respectively). In
conclusion we have the following proposition.
\begin{proposition}\label{prop-32}
Let $A \in \mathcal{A}_d$ such that $P(A)=\Theta_r$, where
$P(z)=z^r-a_0z^{r-1}-\dots-a_{r-1}$. Then, the unique solution of
the matrix differential equation \eqref{Equadiff-S3}, with the
prescribed initial values $X(0)$ and $X'(0)$, is given by
$X(t)=C(t)X(0)+S(t)X'(0)$ with
\begin{equation}\label{solu-S33}
C(t)=\sum_{k=0}^{r-1}\phi_{2k+1}(t)A^k, \quad
S(t)=\sum_{k=0}^{r-1}\phi_{2k}(t)A^k,
\end{equation}
where the $\phi_k(t)$ ($0\leq k\leq r-1$) are described by
\eqref{solu-S32}.
\end{proposition}
Proposition \ref{prop-32} represents a new formulation of the result
of Kolodner (see \cite{ko}). As an application of the last
propositions, let us consider the following example.
\begin{example}\label{example-1} \rm
Let $A\in \mathcal{A}_d$ satisfying
$P(A)=A^2 -1_d=\Theta_d$ (or equivalently $A^2=1_d$),
then $\rho_0(2n)=1, \rho_0(2n+1)=0$ and $\rho_1(2n)=0, \rho_1(2n+1)=0$
for $n\geq 2$.
And a straightforward verification, using Propositions \ref{prop-31}
or \ref{prop-32}, shows that we have
$X(t)=\left(C_0(t)X(0)+S_0(t)X'(0)\right)1_d$ +
$\left(C_1(t)X(0)+S_1(t)X'(0)\right)A$, where
\begin{gather*}
C_0(t)=\frac{1}{2}{\left(ch(t)+\cos(t)\right)} ,\quad
S_0(t)=t+ \frac{1}{2}{\left(sh(t)+\sin(t)\right)},\\
C_1(t)=\frac{1}{2}{\left(ch(t)-\cos(t)\right)} ,\quad
S_1(t)=t+ \frac{1}{2}{\left(sh(t)-\sin(t)\right)}.
\end{gather*}
\end{example}
\subsection{Higher order linear differential equations of Apostol type}
In this Subsection we extend results of Subsection
\ref{Subsection-3.1} to the class of {\it Apostol} linear matrix
differential equations of order $p\geq 2$,
\begin{equation}\label{Equadiff-p}
X^{(p)}(t)=AX(t),
\end{equation}
whose prescribed initial conditions are $X(0)$, $X'(0)$, $\dots$,
$X^{(p-1)}(0)$, where $A \in \mathcal{A}_d$ satisfying the
polynomial equation $P(A)=\Theta_d$, where
$P(z)=z^r-a_0z^{r-1}-\dots-a_{r-1}1_d$. Consider the vector column
$Z(t)=^t(X^{(p-1)}(t), \dots, X(t))$ and the square matrix
$dp\times dp$,
\begin{equation}\label{matrix-B1}
B= \begin{pmatrix}
\Theta_d & \Theta_d & \dots & \Theta_d & A\\
1_d & \Theta_d & \dots & &\Theta_d \\
\Theta_d & 1_d &\Theta_d &\dots &\Theta_d\\
\vdots &\ddots &\ddots &\ddots &\vdots\\
\Theta_d &\dots &\Theta_d &1_d &\Theta_d
\end{pmatrix}.
\end{equation}
It is well known that the solution of the linear matrix differential
equation \eqref{Equadiff-p} can be derived from the usual solution,
\begin{equation}\label{solu-matrix}
Z(t)=e^{tB}Z(0)\quad \text{where }
Z(0)=^t(X^{(p-1)}(0),\dots,X'(0), X(0)),
\end{equation}
of the first order linear system (\ref{Equadiff-12}) defined by
the matrix (\ref{matrix-B1}). That is, we focus our task on the
computation of ${e^{tB}}$, which is essentially
based on the explicit formula of the powers $B^n$($n\geq 0$).
Indeed, by induction we verify that, for every $n\geq 0$, we have
\begin{equation}\label{matrix-B2}
B^{pn}= \begin{pmatrix}
A^n &\Theta_d &\dots & \Theta_d\\
\Theta_d &\ddots &\ddots &\vdots\\
\vdots &\ddots &\ddots &\Theta_d\\
\Theta_d &\dots& \Theta_d &A^n
\end{pmatrix},\quad
B^{pn+i}= \begin{pmatrix}
\Theta_d & \dots & A^{n+1} &\dots & \Theta_d\\ \vdots
&\ddots &
\ddots &\ddots &\vdots\\ \Theta_d & \Theta_d &\dots &\Theta_d &A^{n+1}\\
A^n &\Theta_d &\ddots &\ddots &\vdots\\ \vdots &\ddots
&\ddots &\ddots &\vdots\\
\Theta_d & \dots &A^n &\dots &\Theta_d
\end{pmatrix},
\end{equation}
for $1\leq i\leq p-1$, where the power $A^{n+1}$ of the first line
is located at the $(p-i+1)^{th}$ column and the power $A^n$ of the
first column is located at the $(i+1)^{th}$ line. That is, the
matrix $B^{np+i}=(E_{i,j})_{1\leq i,j\leq p}$ such that $1\leq i\leq
p-1$ and $E_{i,j}\in \mathcal{A}_d$, is described as follows :
$E_{j,p-i+j}=A^{n+1}$ for $j=1,\dots ,i+1$, $E_{i+j,j}=A^{n+1}$ for
$j=1,\dots ,p-i-1$ and $E_{i,j}=\Theta_d$ otherwise. The solution
$X(t)$ of the matrix differential equation \eqref{Equadiff-p}, is
derived by a direct computation. More precisely, the natural
generalization of Proposition \ref{prop-31} is formulated as
follows.
\begin{theorem}\label{thm-34}
Let $A$ be in $\mathcal{A}_d$. Then, the linear matrix
differential equation \eqref{Equadiff-p} of order $p\geq 2$, with
the prescribed initial conditions $X(0)$, $X'(0)$, $\dots$,
$X^{(p-1)}(0)$ has a unique solution $X(t)$ satisfying
$X(t)=C_0(t,A)X(0)+C_1(t,A)X'(0)+\dots +C_{p-1}(t, A)X^{(p-1)}(0)$,
where
\begin{equation}\label{solu-Cj-p}
C_j(t,A)=\sum_{n=0}^{+\infty}\frac{t^{pn+j}}{(np+j)!}A^n,\quad
\text{for every } j \, (0\leq j\leq p-1).
\end{equation}
\end{theorem}
Expression (\ref{matrix-B2}) of the powers of the matrix
(\ref{matrix-B1}), can be easily derived using the general method
of Section \ref{section4}, where this later is based on the
extension of the technique of \cite{br4, mrw}. On the other hand,
properties of the family of the functions $C_j(t,A)$ ($0\leq j\leq
p-1$), may be obtained from (\ref{solu-Cj-p}). Indeed, a simple
verification gives the corollary.
\begin{corollary}\label{coro-35}
Under the hypothesis of Theorem \ref{thm-34}, the functions $C_j(t)$
($0\leq j\leq p-1$) defined by (\ref{solu-Cj-p}) satisfy the
following differential relations,
$$
D^kC_j(t,A)=C_{j-k}(t,A), \; \text{for }\; 0\leq k\leq j\leq p-1,\;
\text{and }\;D^kC_j(0)=\delta_{j,k},
$$
where $D=d/dt$. Moreover, for every $j$ $(0\leq j\leq p-1)$ we have
$C_j(t)=D^{p-j-1}C_{p-1}(t)$.
\end{corollary}
A direct application of (\ref{powers-1})-\eqref{rho-2} and Theorem
\ref{thm-34}, permits us to establish that solutions of the linear
matrix differential equation \eqref{Equadiff-p} are expressed in
terms of the coefficients $a_j$ ($0\leq j\leq r-1$) of the
polynomial $P(z)=z^r-a_0z^{r-1}-\dots-a_{r-1}$ satisfying
$P(A)=\Theta_d$. More precisely, replacing $A^n$ in
(\ref{solu-Cj-p}) by its expression given in (\ref{powers-1})
yields the following proposition.
\begin{proposition}\label{prop-36}
Let $A$ be in $\mathcal{A}_d$ such that $P(A)=\theta_d$, where
$P(z)=z^r-a_0z^{r-1}-\dots-a_{r-1}$. Then, the solution of the
linear differential equation \eqref{Equadiff-p}, with the prescribed
initial conditions $X(0)$, $X'(0)$, $\dots$, $X^{(p-1)}(0)$, is
\[
X(t)=\sum_{k=0}^{r-1}\big[
\sum_{j=0}^{p-1}C_{j,k}(t)X^{(j)}(0)\big]A^k
\]
such that
\[
C_{j,k}(t)=\frac{t^{pk+j}}{(pk+j)!}+\sum_{n=r}^{+\infty}
\frac{t^{pn+j}}{(pn+j)!}\rho_k(n),
\]
with $\rho_k(n)=\sum_{j=0}^ka_{r-k+j-1}\rho(n-j,r)$, where
$\rho(n,r)$ is given by \eqref{rho-2}.
\end{proposition}
Now we obtain an explicit formula for the solution of the
linear matrix differential equation \eqref{Equadiff-p}, in terms
of the roots of the polynomial
$P(z)=z^r-a_0z^{r-1}-\dots-a_{r-1}$ satisfying $P(A)=\Theta_d$.
To this aim, in Expression (\ref{solu-Cj-p}) we substitute for
$A^n$ ($n\geq r$) its expression set forth in Proposition
\ref{prop-23}. Then, a hard straightforward computation leads us
to the following result.
\begin{theorem}\label{thm-37}
Let $A$ be in $\mathcal{A}_d$ such that
$P(A)=\prod_{j=1}^s(A-\lambda_j1_d)^{m_j}=\Theta_d$
($\lambda_i\neq \lambda_j$ for $i\neq j$). Then, the unique
solution of the linear
matrix differential equation \eqref{Equadiff-p}, with the prescribed
initial conditions $X(0)$, $X'(0)$, $\dots $, $X^{(p-1)}(0)$, is
$ X(t)= \sum_{i=0}^{p-1}C_{i}(t,A)X^{(i)}(0)$ such
that
$$
C_{i}(t,A)=\sum_{j=1}^s\varphi_j(t; A),\text{with
} \varphi_j(t; A)=\sum_{k=0}^{m_j-1}
\frac{d^kV_i}{dz^k}(t,\lambda_j)\rho_{jk}(A),
$$
where $V_i(t,z)=\sum_{n=0}^{+\infty}\frac{t^{pn+i}}{(pn+i)!}z^n$ and
$$
\rho_{j,k}(z)=\frac{1}{k!}\Big(\sum_{d=0}^{m_j-k-1}
\alpha_{d,j}(z-\lambda_j)^{k+d}\Big)
\Big(\prod_{t=1, t\neq j}^s\frac{(z-\lambda_t)^{m_t}}
{(\lambda_j-\lambda_t)^{m_t}}\Big),
$$
with $\alpha_{d,j}$ given by \eqref{alpha-3}.
\end{theorem}
The expression of the $C_j(t,A)$ ($0\leq j\leq s$) given in
Theorem \ref{thm-37} seems quite complicated, yet its application
is a powerful tool in many important practical situations as shown
is the following corollaries and Example \ref{example-2}, it also
leads to have explicit formulas. Indeed, when
$p(A)=\prod_{j=1}^r(A-\lambda_j1_d)=\Theta_d$ ($\lambda_i\neq
\lambda_j$, then we show that $\rho_{j0}=\alpha_{0j}=1$, and the
following corollary is derived.
\begin{corollary}\label{coro-38}
Let $A$ be in $\mathcal{A}_d$ such that
$p(A)=\prod_{j=1}^r(A-\lambda_j1_d)=\Theta_d$ ($\lambda_i\neq
\lambda_j$ for $i\neq j$). Then, the unique solution of Equation
\eqref{Equadiff-p}, with the prescribed initial conditions $X(0)$,
$X'(0)$, $\dots$, $X^{(p-1)}(0)$ is $
X(t)=\sum_{i=0}^{p-1}C_{i}(t,A)X^{(i)}(0)$, where
\[
C_i(t,z)=\sum_{j=1}^r \Big(\sum_{n=0}^{+\infty}
\frac{t^{pn+i}}{(pn+i)!}\lambda_j^n\Big)\prod_{t=1,
t\neq j}^r\frac{(z-\lambda_t)}{(\lambda_j-\lambda_t)}.
\]
\end{corollary}
In the particular case where $p=2$, the result of Apostol (see
\cite{ap}) is derived as a simple consequence of the above
corollary. Another important result of Theorem \ref{thm-37} can be
established when $p(A)=(A-\lambda_11_d)^r=\Theta_d$. Indeed, as it
was noticed above, in this case we have $\alpha_{01}=1$ and
$\alpha_{ij}=0$ for each $i\neq 0$ , whence $
\rho_{1,k}(z)=\frac{1}{k!}(z-\lambda_1)^{k}$.
\begin{corollary}\label{coro-39}
Suppose that $A\in\mathcal{A}_d$ satisfies
$p(A)=(A-\lambda_11_d)^r=\Theta_d$. Then, the unique solution of
Equation \eqref{Equadiff-p}, with the prescribed initial conditions
$X(0)$, $X'(0)$, $\dots$, $X^{(p-1)}(0)$, is $ X(t)=
\sum_{i=0}^{p-1}C_{i}(t,A)X^{(i)}(0)$, where
\[
C_{i}(t,A)=\sum_{j=1}^r\frac{1}{(j-1)!}
\frac{d^{j-1}V_i}{d\lambda^{j-1}}(t,\lambda_1)(A-\lambda_11_d)^{j-1}.
\]
\end{corollary}
\begin{example}\label{example-2} \rm
Let $A\in \mathcal{A}_d$ such that $P(A)=\Theta_d$, where
$P(z)=(z-\lambda)^2(z-\mu)$ (with $\lambda \neq \mu$). Then, the
solution of the third-order matrix differential equation
$X'''(t)=AX(t)$, with the prescribed initial conditions $X(0)$,
$X'(0)$ and $X''(0)$, is given by
$X(t)=C_0(t,A)X(0)+C_1(t,A)X'(0)+C_2(t,A)X''(0)$. And by Corollary
\ref{coro-38} we have $C_1(t,A)=DC_2(t,A)$, $C_0(t,A)$
=$D^2C_2(t,A)$ ($D=d/dt$), where
$$
C_2(t,A)=V_2(t,\lambda)\big[1_d+\frac{A-\lambda1_d
}{\mu-\lambda}\big]\phi_1(A)+\frac{dV_2}{dz}(t,\lambda)(A-\lambda1_d
)\phi_1(t, \lambda)+V_2(t,\mu)\phi_2(A),
$$
with
\[
V_2(t,z)=\sum_{n=0}^{+\infty}\frac{t^{3n+2}}{(3n+2)!}z^n,\quad
\phi_1(z)=\frac{z-\mu}{\lambda-\mu}
\]
and
${\phi_2(z)=\frac{\left(z-\lambda\right)^2}{\left(\mu-\lambda\right)^2}}$.
\end{example}
\begin{remark} \label{rmk3.11}
As shown before, our method for studying these kinds of matrix
differential equations is more direct and does not appeal to other
technics. Meanwhile, it seems for us that the method of Verde Star
based on the technics of divided differences can be also applied
for studying such equations (see \cite{vs}).
\end{remark}
\section{Solutions of Equation \eqref{Equadiff-1}:
combinatorial setting}\label{section4}
We are interested here in the combinatorial solutions of \eqref{Equadiff-1}, where the
exponential generating function of the family of sequences
$\{\rho(n-j,r)\}_{n\geq 0}$ ($0\leq j\leq r-1$) is defined by
\eqref{rho-1}.
Let $A_0,\dots, A_{r-1}$ and $S_0,\dots, S_{r-1}$ be two finite
sequences of $\mathcal{A}_d$ such that $A_{r-1}\neq 0$. Let
${\mathcal{ C}}^{\infty}(\mathbb{R}\ ; \mathcal{A}_d)$ be the
$\:\mathbb{C}$-vector space of functions of class
${\mathcal{C}}^{\infty}$. Consider the class of linear matrix
differential equations of higher-order \eqref{Equadiff-1}, with solutions in
$X\in {\mathcal{ C}}^{\infty}(\mathbb{R}\ ; \mathcal{A}_d)$ and subject to the
initial conditions, $X(0)$,
$X'(0)$, $\dots $, $X^{(r-1)}(0)$. Set $Z(t)=^t(X^{(r-1)}(t),
\dots, X(t))$ and $Z(0)=^t(X^{(r-1)}(0), \dots, X(0))$. A
standard computation shows that \eqref{Equadiff-1} is equivalent
to the following first-order matrix differential equation,
\begin{equation}\label{comp-equadiff-1}
Z'(t)=BZ(t),
\end{equation}
where $B \in \mathcal{M}(r,\;\mathcal{A}_d)$ is the
companion matrix
\begin{equation}\label{comp-matr}
B= \begin{pmatrix} A_0 & A_1 & \dots & & A_{r-1}\\
1_d & \Theta_d & \dots & &\Theta_d \\
\Theta_d & 1_d &\Theta_d &\dots &\Theta_d\\
\vdots &\ddots &\ddots &\ddots &\vdots\\
\Theta_d &\dots &\Theta_d &1_d &\Theta_d
\end{pmatrix}.
\end{equation}
It is well known that the solution of the linear matrix differential
equation \eqref{Equadiff-1} is
$ Z(t)=e^{tB}Z(0)$, where the expression of $ e^{tB}$ depends on the
computation of the powers $B^n$ with the aid of $A_0,\dots, A_{r-1}$.
To this aim, we apply the recent technique for computing
the powers of the usual companion matrix (\ref{comp-matr}) (see
\cite{br4, mrw}). Indeed, let $\{Y_{n,s}\}_{n\geq 0}$ ($0\leq s\leq
r-1$) be the class of sequences \eqref{recursive-1} in
$\mathcal{A}_d$ such that $ Y_{n,s}= \delta_{n,s}1_d$
($\delta_{n,s}= 1$ if $n=s$ and 0 if not) for $0\leq n\leq r-1$ and
\begin{equation}\label{recursive-s}
Y_{n+1,s}= A_0Y_{n,s}+A_1
Y_{n-1,s}+\dots+A_{r-1} Y_{n-r+1,s}, \quad \text{for } n\geq r-1.
\end{equation}
\begin{lemma}\label{lemm-41}
Let $A_0,\dots, A_{r-1}$ be in $\mathcal{A}_d$ with $A_{r-1}\neq 0$.
Then
\begin{equation}\label{comp-power}
B^n= \begin{pmatrix}
Y_{n+r-1,r-1} & Y_{n+r-1,r-2} & \dots & Y_{n+r-1, 0} \\
Y_{n+r-2,r-1} & Y_{n+r-2,r-2} & \dots & Y_{n+r-2,0} \\
\vdots &\vdots
&\ddots &\vdots\\
Y_{n,r-1} & Y_{n,r-2} & \dots & Y_{n,0} \\
\end{pmatrix},\quad \text{for every } n\geq 0.
\end{equation}
\end{lemma}
Expression (\ref{comp-power}) implies that
$e^{tB}=(C_{i,j}(t))_{0\leq i,\; j\leq r-1}$, where $C_{i,j}(t)$
($0\leq i,\; j\leq r-1$) is the exponential generating functions of
the sequence $\{Y_{n+r-1-i,r-1-j}\}_{n\geq 0}$ ($0\leq i,\; j\leq
r-1$). Thus, the characterization of solutions of the linear matrix
differential equation \eqref{Equadiff-1}, can be formulated as
follows.
\begin{theorem}\label{thm-42}
Let $A_0,\dots, A_{r-1}$ be in $\mathcal{A}_d$ such that
$A_{r-1}\neq 0$. Then, the solution of the linear matrix
differential equation \eqref{Equadiff-1} is
$X(t)=\sum_{j=0}^{r-1}C_j(t)X^{(j)}(0)$, where
$C_j(t)=\sum_{n=0}^{+\infty}Y_{n,j}t^n/n!$ ($0\leq j\leq r-1$) is
the exponential generating function of the sequences
$\{Y_{n,j}\}_{n\geq 0}$ ($0\leq j\leq r-1$) defined by
\eqref{recursive-s}.
\end{theorem}
A result similar to Theorem \ref{thm-42} has been given by Verde
Star under another form, by using the divided difference method
for solving linear differential equations \eqref{Equadiff-1} (see
Section 5 of \cite{vs2}). Indeed, Equation \eqref{Equadiff-1} is
analogous to the Equation (5.6) given by Verde Star whose
solutions (5.13), submitted to the initial conditions
$C_0,C_1,\dots , C_s$, are expressed in terms of a sequence
$\{P_{k+j}\}_{{0\leq j\leq s},k\geq 0}$ (see \cite{vs2}),
satisfying a linear recursive relation analogous to
\eqref{recursive-s}, moreover $\{P_{k+j}\}_{{0\leq j\leq s},k\geq
0}$ is nothing else but the sequence $\{Y_{k,s-j}\}_j$ defined
previously by \eqref{recursive-s}. Comparison of our procedure
with the Verde Star's method leads to infer that the functions
$C_0(t), C_1(t),\dots, C_{r-1}(t)$ of Theorem \ref{thm-42} are
identical to the functions $G_{0,1}(t), G_{0,2}(t),\dots,
G_{0,s}(t)$ considered in \cite{vs2}.
Furthermore when the condition $A_jA_k=A_kA_j$ ($0\leq j,k\leq
r-1$) is satisfied, expressions \eqref{combinat-1}--\eqref{rho-1}
show that the combinatoric formula of sequences
\eqref{recursive-s} can be written explicitly as follows,
\begin{equation}\label{combinat-2}
Y_{n,j}=\sum_{k=0}^jA_{r-j+k-1}\rho(n-k,r),\quad\text{for } n\geq
r,
\end{equation}
where the $\rho(n,r)$ are defined in \eqref{rho-1} for $n\geq r$,
with $\rho(r,r)=1_d$ and $\rho(n,r)=\Theta_d$ for $n\leq r-1$.
Therefore, Expression (\ref{combinat-2}) implies that
$C_j(t)=\frac{t^j}{j!}1_d+\sum_{k=0}^jA_{r-j+k-1}g_k(t)$, where
$g_k(t)=\sum_{n=0}^{+\infty}\rho(n-k,r)t^n/n!$ is the exponential
generating functions of the sequence $\{\rho(n-k,r)\}_{n\geq 0}$.
Moreover, we verify easily that the functions $g_j(t)$
($0\leq k\leq j$) satisfy
${g_j^{(j-k)}(t)=\frac{d^{j-k}g_j}{dt^{j-k}}(t)}
$=$g_k(t)$. Hence, an extension of Proposition \ref{prop-31} and
Theorem \ref{thm-34} can be formulated as follows.
\begin{proposition}\label{prop-43}
Let $A_0,\dots, A_{r-1}$ be in $\mathcal{A}_d$ such that
$A_iA_k=A_kA_i$ for $0\leq i$, $k\leq r-1$.
Then the solution of the
linear matrix differential equation \eqref{Equadiff-1}, subject to
the prescribed initial values $X(0), X'(0), \dots ,$
$X^{(r-1)}(0)$, is the following
\begin{equation}\label{solu-g2}
X(t)=\sum_{j=0}^{r-1}\Delta_j(t)X^{(j)}(0)=
\sum_{j=0}^{r-1}\Big(\frac{t^j}{j!}+\Pi_j(D)g_j(t)\Big)X^{(j)}(0),
\end{equation}
such that
\[
\Delta_j(t)=\frac{t^j}{j!}+\sum_{k=0}^jA_{r-j+k-1}g_k(t),\quad
\Pi_j(D)=\frac{t^j}{j!}+\sum_{k=0}^jA_{r-j+k-1}D^{j-k},
\]
where
$D=d/dt$ and $ g_j(t)$ is the exponential generating function of the
sequence $\{\rho(n-j,r)\}_{n\geq 0}$.
\end{proposition}
Proposition \ref{prop-43} shows that solutions of the linear
matrix differential equation \eqref{Equadiff-1} may be given in terms
of the exponential generating function of the class of
sequences $\{Y_{n,j}\}_{n\geq 0}$ ($0\leq j\leq r-1$) defined by
\eqref{recursive-s}. In the same way, expressions
\eqref{combinat-1}-\eqref{rho-1} may be used to obtain the
combinatoric formulas of the sequences in \eqref{recursive-s}. Moreover,
solutions of the linear matrix differential equation
\eqref{Equadiff-1}, subject to the prescribed initial values
$X(0), X'(0), \dots ,$ $X^{(r-1)}(0)$, can be expressed in terms
the exponential generating functions of the class of sequences
$\{\rho(n-k,r)\}_{n\geq 0}$ ($0\leq k\leq r-1$).
\begin{remark}\label{remark-31}\rm
Consider a unitary topological
$\mathbb{C}$-algebra ${\mathcal{A}}$ instead of the $\mathbb{C}$-algebra of the
square matrices $\mathcal{A}_d$. Suppose that $A\in
{\mathcal{A}}$ is an algebraic element satisfying $P(A)=0$, where
$P(z)=z^r-a_0z^{r-1}-\dots-a_{r-1}$ ($a_j\in \mathbb{C}$ for $0\leq j\leq
r-1$). Then, all results of the preceding Sections are still
valid. On the other hand, it's easy to show that Theorem
\ref{thm-34} and Corollary \ref{coro-35}, on the solutions of the
linear matrix differential equation \eqref{Equadiff-p}, do not
depend on the condition that $A$ is algebraic.
\end{remark}
\begin{remark}[Future problem] \label{rmk4.5} \rm
Solutions of the linear matrix differential equation
\eqref{Equadiff-3}, can be also described using some recursive
relations and the exponential generating functions of the
combinatorial sequences \eqref{recursive-1}. One of the main
problems is to study the spectral aspect of solutions of these
classes of differential equations. More precisely, suppose that
the two (non trivial) matrices $A_0$, $A_1$ appearing in
\eqref{Equadiff-3} satisfy the following polynomial equations
$P_1(A_0)=P_2(A_1)=\Theta_d$, where
$P_1(z)=\prod_{j=0}^{s_1}(z-\lambda_j)^{p_j}$ and
$P_2(z)=\prod_{j=0}^{s_2}(z-\mu_j)^{q_j}$. The problem can be
formulated as follows : {\it study the solutions of the linear
matrix differential equation \eqref{Equadiff-3}, in terms of the
$\lambda_j$ ($0\leq j\leq s_1$) and $\mu_j$ ($0\leq j\leq s_2$)}.
Some investigations are currently done for this purpose.
\end{remark}
\subsection*{Acknowledgments}
The authors would like to express
their sincere gratitude to Professor L. Verde Star for his helpful
suggestions and fruitful correspondence that improved this paper.
They are also grateful to the anonymous referee for his/her valuable
comments on an earlier version of this paper.
\begin{thebibliography}{99}
\bibitem{ap} T.M. Apostol;
\emph{Explicit formulas for solutions of the second order matrix differential equation
$Y''=AY$,} Amer. Math. Monthly 82 (1975), pp. 159-162.
\bibitem{br1} R. Ben Taher and M. Rachidi;
\emph{Linear recurrence relations in the algebra of matrices and applications,}
Linear Algebra and Its Applications, Vol. 330 (1-3) (2001), pp. 15-24.
\bibitem{br2} R. Ben Taher and M. Rachidi;
\emph{Some explicit formulas for the polynomial decomposition of the matrix
exponential and applications,} Linear Algebra and Its Applications,
350, Issue 1-3 (2002), pp. 171-184.
\bibitem{br3} R. Ben Taher, I. Bensaoud and M. Rachidi;
\emph{Dynamic solutions and computation of the powers and
exponential of matrix,} International Journal of Mathematics, Game
Theory, and Algebra 13 (2003), No. 6, p. 455-463.
\bibitem{br4} R. Ben Taher and M. Rachidi;
\emph{On the matrix powers and exponential by the $r$-generalized
Fibonacci sequences methods. The companion matrix case}, Linear
Algebra and Applications, 370 (2003), pp. 341-353.
\bibitem{cy} H.W. Cheng and S.S.-T. Yau;
\emph{On more explicit formulas for matrix exponential,}
Linear Algebra Appl. 262 (1997), pp. 131-163.
\bibitem{mrw} B. El Wahbi, M. Mouline and M. Rachidi;
\emph{Solving nonhomogeneous recurrence relations
by matrix methods,} The Fibonacci Quarterly, 40-2 (2002), pp.
109-117.
\bibitem{ga} F. R. Gantmacher;
\emph{Theory of matrices}, Chelsea Publishing Company, New York (1959).
\bibitem{ko} I. I. Kolodner;
\emph{On $exp(tA)$ with $A$ satisfying a polynomial},
J. Math.Anal. and Appl. 52 (1975), pp. 514-524.
\bibitem{le} C. Levesque;
\emph{On the $m-th$ order linear recurrences}, Fibonacci Quarterly 23 (4),
(1985) : 290-295.
\bibitem{leo} I. E. Leonardo;
\emph{The matrix exponential}, SIAM Review Vol. 38, No. 3 (1996),
pp. 507-512.
\bibitem{mv} C. Moler and C. Van Loan;
\emph{Nineteen dubious ways to compute the exponential of matrix,
twenty-Five years later,} SIAM Review Vol.45 (1), (2003), pp.
3-49.
\bibitem{mr1} M. Mouline and M. Rachidi;
\emph{Application of Markov Chains properties to $r$-Generalized
Fibonacci Sequences,} Fibonacci Quarterly 37 (1999), pp. 34-38.
\bibitem{mr2} M. Mouline and M. Rachidi;
\emph{Suites de Fibonacci g\'en\'eralis\'ees, Th\'eor\`eme de Cayley-Hamilton
et chaines de Markov,} Rendiconti del Seminario Matematico di
Messina Serie II, Tomo XIX, No. 4 (1996/97), pp. 107-115.
\bibitem{pu} E. J. Putzer;
\emph{Avoiding the Jordan canonical form in the discussion of linear systems with constant
coefficients,} Amer. Math. Monthly 73 (1966), pp. 2-7.
\bibitem{rp} Hans-J. Runckel and U. Pitelkow;
\emph{Practical Computation of matrix functions,} Linear Algebra Appl.
49 (1983), pp. 161-178.
\bibitem{vs} L. Verde-Star;
\emph{Operator identities and the solution of linear matrix difference
and differential equations,} Studies in Applied Mathematics 91 (1994), pp.
153-177.
\bibitem{vs1} L. Verde-Star;
\emph{Divided differences and linearly recurrent sequences,}
Stud. Appl. Math. 95 (1995), pp. 433-456.
\bibitem{vs2} L. Verde-Star;
\emph{Solutions of linear differential equations by the
method of divided differences,} Adv. in Applied Math. 16 (1995), pp.
484-508.
\end{thebibliography}
\end{document}