\documentclass[reqno]{amsart} \AtBeginDocument{{\noindent\small {\em Electronic Journal of Differential Equations}, Vol. 2003(2003), No. 79, pp. 1--18.\newline ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu \newline ftp ejde.math.swt.edu (login: ftp)} \thanks{\copyright 2003 Southwest Texas State University.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE--2003/79\hfil Differential transfer matrix method] {Analytical solution of linear ordinary differential equations by differential transfer matrix method} \author[Sina Khorasani \& Ali Adibi \hfil EJDE--2003/79\hfilneg] {Sina Khorasani \& Ali Adibi} \address{School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA 30332-0250, USA} % Fax: (404) 894-4641} \email[Sina Khorasani]{sina.khorasani@ece.gatech.edu} \date{} \thanks{Submitted June 2, 2003. Published August 5, 2003.} \subjclass[2000]{34A05, 34A25, 34A30} \keywords{Ordinary differential equations, differential transfer matrix method, \hfill\break\indent matrix exponential} \begin{abstract} We report a new analytical method for finding the exact solution of homogeneous linear ordinary differential equations with arbitrary order and variable coefficients. The method is based on the definition of jump transfer matrices and their extension into limiting differential form. The approach reduces the $n$th-order differential equation to a system of $n$ linear differential equations with unity order. The full analytical solution is then found by the perturbation technique. The important feature of the presented method is that it deals with the evolution of independent solutions, rather than its derivatives. We prove the validity of method by direct substitution of the solution in the original differential equation. We discuss the general properties of differential transfer matrices and present several analytical examples, showing the applicability of the method. \end{abstract} \maketitle \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \numberwithin{equation}{section} \allowdisplaybreaks \section{Introduction} Solution of ordinary differential equations (ODEs) is in general possible by different methods [1]. Lie group analysis [2-6] can be effectively used in cases when a symmetry algebra exists. Factorization methods are reported for reduction of ODEs into linear autonomous forms [7,8] with constant coefficients, which can be readily solved. But these methods require some symbolic transformations and are complicated for general higher order ODEs, and thus not suitable for many computational purposes where the form of ODE is not available in a simple closed form. General second-order linear ODEs can be solved in terms of special function [9] and analytic continuation [10]. Spectral methods [11,12] and orthogonal polynomial [13] solutions of linear ODEs have been also reported. Solution properties of linear ODEs with arbitrary orders are discussed in a report [14], and several approximation methods have been also published [15]. All of these methods express the solution in terms of infinite series, relying mainly on numerical computations. Analytic solution of such equations has been also recently presented [16], but only at a singular point. Linear ODEs of the $n$th-order can also be transformed to a system of $n$ linear first-order differential equations by reduction method [17], where the derivative of one parameter is regarded as the next parameter. This approach is based on the evolution of the solution function and its derivatives up to the order $n$. This method is therefore not useful when the evolution of independent solutions is of interest rather than the derivatives, such as waves in inhomogeneous media. Recently, we have developed a new matrix-based exact analytical method in order to study the propagation of electromagnetic waves in non-uniform media with an arbitrary refractive index function [18,19]. The method is based on the modification of the well-known transfer matrix theory in optics [20] and quantum mechanics [21], and its extension into differential form. This method is simple, exact, and efficient, while reflecting the basic properties of the physical system described by the differential equation. We have also developed proper formalism for anisotropic non-uniform media, in which the governing equations are no longer scalar [22]. In this method, the solution is found by integration of an exponent square matrix and then taking its matrix exponential. In this paper, we show that the differential transfer matrix method (DTMM) [18,19] can be also used for exact solution of linear ODEs with arbitrary variable coefficients. We only consider the homogeneous equation and its linear independent solutions, since once the linear independent solutions are known, the particular solution can be found by the method of variation of parameters due to Lagrange [2]. We show that the $n$th order ODE is transformed into a system of $n$ homogeneous first order linear equations, which can be integrated either numerically by Runge-Kutta [23,24] or reflexive [25] methods, or analytically by matrix exponential techniques. Through direct substitution we rigorously show that the presented analytical solution satisfies the corresponding differential equation, and justify the approach through several illustrative examples. We also show that the famous Abel-Liouville-Ostogradsky theorem for the Wronskian of a linear ODE [26-28] can be easily proved through the DTMM. In $\S 2$ we describe the outline of the problem and then formulate the jump transfer matrix method in $\S 3$. In $\S 4$ differential formulation, the fundamental theorem of DTMM, and treatment of singularities are described. Finally, the application of the DTMM to several ODEs of various orders is presented. \section{Statement of Problem} Consider the homogeneous ordinary linear differential equation of order $n$ $$\label{eq1} \mathbb{L} f(x) = 0, \quad f:\mathbb{S} \mapsto \mathbb{C}$$ where $f(x)$ is an unknown analytic complex function in the set of complex variables $\mathbb{C}$ with the connected domain $\mathbb{S} \subset \mathbb{C}$, and $\mathbb{L}$ is a linear operator given by $$\label{eq2} \mathbb{L} = \sum_{m = 0}^n {a_m (x)\frac{d^m}{dx^m}} , \quad a_n (x) \equiv 1.$$ Here we assume that $a_m :\mathbb{S} \mapsto \mathbb{C},m = 0,\ldots, n - 1$ are arbitrary analytic complex functions of $x$, being referred to as variable coefficients. If coefficients are constant, then under the condition described just below, a general solution to (\ref{eq1}) is [1,29] $$\label{eq3} f(x) = c_1 \exp \left( {k_1 x} \right) + c_2 \exp \left( {k_2 x} \right) + \cdots + c_n \exp \left( {k_n x} \right),$$ where $c_i ,i = 1,\ldots, n$ are complex constants determined from either boundary or initial conditions, and $k_i ,i = 1,\ldots, n$, being referred here to as wavenumbers, satisfy the characteristic equation $$\label{eq4} \sum_{m = 0}^n {a_m k_i ^m} = 0,\quad i = 1,\ldots, n.$$ Moreover, for (\ref{eq3}) to be a general solution of (\ref{eq1}) with constant coefficients, it is necessary that $\forall i,j;k_i \ne k_j$, i.e. (\ref{eq4}) has exactly $n$ distinct roots. Now, if coefficient functions are not constant, then we propose a general solution to (\ref{eq1}) of the form $$\label{eq5} f(x) = f_1 (x)\exp \left[ {k_1 (x) x} \right] + f_2 (x) \exp \left[ {k_2 (x) x} \right] + \cdots + f_n (x)\exp \left[ {k_n (x)x} \right],$$ or equivalently \begin{gather} \label{eq6} f(x) = \exp \left[ {{\rm {\bf \Phi }}(x)} \right]^{\,t}\,{\rm {\bf F}}(x),\\ \label{eq7} {\rm {\bf F}}(x) = \left[ \begin{matrix} f_1 (x) \\ f_2 (x) \\ \vdots \\ f_n (x) \end{matrix} \right],\quad {\rm {\bf \Phi }}(x) = \left[ \begin{matrix} k_1 (x)x \\ k_2 (x)x \\ \vdots \\ k_n (x)x \end{matrix}\right], \end{gather} where the superscript $t$ denotes the transposed matrix. Hereinafter, we shall refer to ${\rm {\bf F}}(x)$ and ${\rm {\bf \Phi }}(x)$ respectively as envelope and phase vector functions. In the above equation $f_i (x)$ are functions to be determined and $k_i(x)$, being referred to as wavenumber functions, are complex functions of $x$ satisfying the generalized algebraic characteristic equation $$\label{eq8} \sum_{m = 0}^n {a_m (x)\,k_i ^{\,m}(x)} = 0,\quad i = 1,\ldots, n.$$ For the sake of convenience, we here define the exponential of a vector ${\rm {\bf v}} = \{ {v_i }\}_{n\times 1}$ as $\exp ({\rm {\bf v}}) = \{ {\exp ( {v_i } )} \}_{n\times 1}$. In contrast, $\exp ( {\rm {\bf M}} )$ when \textbf{M} is a square matrix represents the matrix exponential of \textbf{M}, given by $$\label{eq8a} \exp \left( {\rm {\bf M}} \right) = {\rm {\bf I}} + \sum_{m = 1}^\infty {\frac{1}{m!}{\rm {\bf M}}^m}.$$ We define a domain $\mathbb{S}$ to be non-degenerate if (\ref{eq8}) has exactly $n$ distinct roots for all $x \in \mathbb{S}$, and $m$-fold degenerate if it is non-degenerate everywhere, but at finite number of isolated points, being referred to as singularities, at which (\ref{eq8}) has a minimum of $n-m$ distinct roots. If $\mathbb{S}$ is at least 1-fold degenerate for all $x \in \mathbb{S}$, then $\mathbb{S}$ is referred to as entirely degenerate. Later, we shall observe that non-singular solutions require $\mathbb{S}$ to be non-degenerate, and degenerate domains require special treatment. We also define the diagonal wavenumber matrix as $$\label{eq9} {\rm {\bf K}}(x) = \left[ {k_i (x)\delta _{ij} } \right]_{\,n\times n} = {\rm diag} \left[ k_1(x), \ldots, k_n(x) \right],$$ with $\delta _{ij}$ being the Kronecker delta. Obviously, the diagonal elements of ${{\rm {\bf K}}(x)}$, and thus ${{\rm {\bf K}}(x)}$ itself must satisfy (\ref{eq8}). Moreover, $\left| {{\rm {\bf K}}(x)} \right| = a_0 (x)$ and $\mathop{\rm tr}\left\{ {{\rm {\bf K}}(x)} \right\} = a_{n - 1} (x)$, where $\mathop{\rm tr}\left\{\cdot\right\}$ denotes the trace operator. The phase vector and wavenumber matrix are thus related as ${\rm {\bf \Phi }}(x) = x{\rm {\bf K}}(x)\,{\rm {\bf 1}}_{n\times 1}$, where ${\bf 1}_{n\times 1}$ is a vector with unity elements. Hence, one has $\exp \left[ {{\rm {\bf \Phi }}(x)} \right] = \exp \left[ {x{\rm {\bf K}}(x)} \right]\,{\rm {\bf 1}}_{n\times 1}$. \section{Jump Transfer Matrices} \subsection{Transfer Matrix of a Single Jump} Without loss of generality, suppose that $\mathbb{S} \subset \mathbb{R}$ with $\mathbb{R}$ being the set of real numbers. Furthermore, we let the variable coefficients $a_i(x),i=1,\ldots, n$ to be such stepwise constant discontinuous functions at $x=X$ that ${\rm {\bf K}}\left( {x < X} \right) = {\rm {\bf K}}_A$ and ${\rm {\bf K}}\left( {x > X} \right) = {\rm {\bf K}}_B$ become constant matrices. We furthermore accept that both the sub-spaces $\mathbb{S}_A = \left\{ { x : x \in \mathbb{S},x < X} \right\}$ and $\mathbb{S}_B = \left\{ { x : x \in \mathbb{S},x > X} \right\}$ are non-degenerate. In this case, the solution in $\mathbb{S}$ would be given by $$\label{eq10} f(x) = \begin{cases} \exp \left[ {{\rm {\bf \Phi }}_A (x)} \right]^{\,t}\,{\rm {\bf F}}_A , & x \in \mathbb{S}_A \\ \exp \left[ {{\rm {\bf \Phi }}_B (x)} \right]^{\,t}\,{\rm {\bf F}}_B , & x \in \mathbb{S}_B \,, \end{cases}$$ where the envelope vectors ${\rm {\bf F}}_A$ and ${\rm {\bf F}}_B$ must be constant, and ${\rm {\bf \Phi }}_A (x) = \left\{ {_A k_i x} \right\}_{n\times 1}$, ${\rm {\bf \Phi }}_B (x) = \left\{ {_B k_i x} \right\}_{n\times 1}$ with ${ }_Ak_i$ and ${ }_Bk_i$ satisfying the characteristic equation (\ref{eq4}) in $\mathbb{S}_A$ and $\mathbb{S}_B$. Analyticity of $f(x)$ across $x = X$ requires that $$\label{eq11} f^{(m)}( {X^ - }) = f^{(m)}( {X^ + }),\quad m = 0,\ldots, n - 1,$$ where $f^{(m)}(x)$ is the $m$th derivative of $f(x)$ with respect to $x$. Expanding (\ref{eq11}) results in the set of equations $$\label{eq12} {\rm {\bf D}}_B \exp \left( {X{\rm {\bf K}}_B } \right){\rm {\bf F}}_B = {\rm {\bf D}}_A \exp \left( {X{\rm {\bf K}}_A } \right){\rm {\bf F}}_A,$$ in which ${\rm {\bf K}}_A$ and ${\rm {\bf K}}_B$ are given in (\ref{eq9}). Also ${\rm {\bf D}}_A = \left[ {{ }_Ak_j ^{i - 1}} \right]_{\,n\times n}$ and ${\rm {\bf D}}_B = \left[ {{ }_Bk_j ^{i - 1}} \right]_{\,n\times n}$ are matrices of Vandermonde type [30] given by $$\label{eq13} {\rm {\bf D}}_A = \left[ \begin{matrix} 1 & 1& \cdots & 1 \\ {_A k_1 } & {_A k_2 } & \cdots & {_A k_n } \\ \vdots & \vdots & & \vdots \\ {_A k_1 ^{n - 1}} & {_A k_2 ^{n - 1}} & \cdots & {_A k_n ^{n - 1}} \end{matrix} \right], \quad {\rm {\bf D}}_B = \left[ \begin{matrix} 1 & 1 & \cdots & 1 \\ {_B k_1 } & {_B k_2 } & \cdots & {_B k_n } \\ \vdots & \vdots & & \vdots \\ {_B k_1 ^{n - 1}} & {_B k_2 ^{n - 1}} & \cdots & {_B k_n ^{n - 1}} \end{matrix} \right].$$ One can rewrite (\ref{eq12}) as \begin{gather} \label{eq14} {\rm {\bf F}}_B = {\rm {\bf Q}}_{A \to B} {\rm {\bf F}}_A,\\ \label{eq15} {\rm {\bf Q}}_{A \to B} = \exp \left( { - X{\rm {\bf K}}_B } \right)\,{\rm {\bf D}}_B ^{ - 1}{\rm {\bf D}}_A \exp \left( {X{\rm {\bf K}}_A } \right), \end{gather} in which ${\rm {\bf Q}}_{A \to B}$ is defined as the jump transfer matrix from $A$ to $B$. Similarly, one can define the downward jump transfer matrix as ${\rm {\bf Q}}_{B \to A} = {\rm {\bf Q}}_{A \to B}^{ - 1}$. The determinant of the jump transfer matrix ${\rm {\bf Q}}_{A \to B}$ is $$\label{eq16} \left| {{\rm {\bf Q}}_{A \to B} } \right| = \exp \left[ {X\left( {\mathop{\rm tr}\left\{ {{\rm {\bf K}}_A } \right\} - \mathop{\rm tr}\left\{ {{\rm {\bf K}}_B } \right\}} \right)} \right]\frac{\left| {{\rm {\bf D}}_A } \right|}{\left| {{\rm {\bf D}}_B } \right|},$$ in which we have used the identity $\left|\exp({\bf M})\right|=\exp(\mathop{\rm tr}\left\{ {\bf M}\right\} )$. The determinants can be also evaluated by noting the fact that $\left| {{\rm {\bf D}}_r } \right|$ are determinants of Vandermonde matrices, given by [30] $$\label{eq17} \left| {{\rm {\bf D}}_r } \right| = \begin{cases} 1, & n = 1 \\ \prod_{i > j} {\left( {{ }_rk_i - { }_rk_j } \right)} , & n > 1 \end{cases} \quad r = A,B.$$ Here, the product operator runs over all possible ordered pairs $(i,j)$ with $i>j$. Now since $\mathop{\rm tr}\left\{ {{\rm {\bf K}}_r } \right\} = \sum_{i = 1}^n {{ }_rk_i } = { }_ra_{n - 1}$ with $r = A,B$, (\ref{eq16}) can be expanded as $$\label{eq18} \left| {{\rm {\bf Q}}_{A \to B} } \right| = \exp \left[ X ({ }_Aa_{n - 1}-{}_Ba_{n - 1})\right]\prod_{i > j} {\frac{\left( {{ }_Ak_i - { }_Ak_j } \right)}{\left( {{ }_Bk_i - { }_Bk_j } \right)}}.$$ \subsection{Transfer Matrix of Multiple Jumps} When variable coefficients $a_i(x),i=1,\ldots ,n$ are stepwise functions of $x$ with arbitrary number of discontinuities or jumps over interfaces, ${\rm {\bf K}}(x)$ can be expressed as $$\label{eq19} {\rm {\bf K}}(x) = {\rm {\bf K}}_p ,\quad x \in \mathbb{S}_p ,\quad p = 0, \ldots, P ,$$ where the subset $\mathbb{S}_p = \left\{ {\left. x \right|X_p < x < X_{p + 1} } \right\},p = 0,\ldots, P$ is referred to as the $p$th layer. Obviously, $X_0 = \inf \left\{\mathbb{S} \right\} = \inf \left\{ {\mathbb{S}_0 } \right\}$ and $X_{P + 1} = \sup \left\{\mathbb{S} \right\} = \sup \left\{ {\mathbb{S}_P } \right\}$. In this case, the corresponding envelope vectors are related as $$\label{eq20} {\rm {\bf F}}_s = {\rm {\bf Q}}_{r \to s} {\rm {\bf F}}_r ,\quad r,s = 0, \ldots, P,$$ where ${\rm {\bf Q}}_{r \to s}$ is the transfer matrix from layer $r$ to $s$, obtained by multiplication of single jump transfer matrices as $$\label{eq21} {\rm {\bf Q}}_{r \to s} = {\rm {\bf Q}}_{s - 1 \to s} {\rm {\bf Q}}_{s - 2 \to s - 1} \cdots {\rm {\bf Q}}_{r \to r + 1}.$$ \subsection{Properties of Jump Transfer Matrix } The transfer matrix ${\rm {\bf Q}}_{r \to s}$ satisfies the basic properties \begin{gather} \label{eq21a} {\rm {\bf Q}}_{r \to r} = {\rm {\bf I}},\quad \mbox{(self-projection)}\\ \label{eq21b}{\rm {\bf Q}}_{r \to s} = {\rm {\bf Q}}_{s \to r}^{ - 1},\quad \mbox{(inversion)}\\ \label{eq21c} {\rm {\bf Q}}_{r \to s} = {\rm {\bf Q}}_{u \to s} {\rm {\bf Q}}_{r \to u},\quad\mbox{(decomposition)}\\ \label{eq21d} \left| {{\rm {\bf Q}}_{r \to s} } \right| = \exp \left[ {\sum_{p = r}^s {X_p \sum_{i = 1}^n {\left( {{ }_pk{ }_i - { }_{p + 1}k{ }_i} \right)} } } \right]\prod_{i > j} {\frac{\left( {{ }_rk_i - { }_rk_j } \right)}{\left( {{ }_sk_i - { }_sk_j } \right)}},\quad\mbox{(determinant)} \end{gather} with $r,s,u = 0, \ldots, P$. The properties (\ref{eq21a}), (\ref{eq21b}), and (\ref{eq21c}) are direct results of the definition of transfer matrices, and the determinant property (\ref{eq21d}) follows (\ref{eq21}) and (\ref{eq18}). From (\ref{eq18}) it can be observed that for systems with vanishing $_r a_{n - 1} ,r = 0, \ldots, P$, the determinant property simplifies into $$\label{eq22} \left| {{\rm {\bf Q}}_{r \to s} } \right| = \prod_{i > j} {\frac{\left( {{ }_rk_j - { }_rk_i } \right)}{\left( {{ }_sk_j - { }_sk_i } \right)}}.$$ This expression has the important feature of being dependent only on the local properties of the initial layer $r$ and target layer $s$. The transfer matrix ${\rm {\bf Q}}_{r \to s}$ satisfies the scaling property, that is it remains unchanged under transformations $X_p \to \alpha X_p$ and ${ }_rk_i \to \alpha ^{ - 1}{ }_rk_i$. This can be directly observed either from (\ref{eq1}) by making the change of variable $x \to \alpha x$, which accordingly scales the roots in an inverse manner, or direct substitution in (\ref{eq15}). It also preserves the shifting property $$\label{eq22a} {\rm {\bf Q}}_{r \to s} = \exp \left( { - \xi {\rm {\bf K}}_s } \right)\;{\rm {\bf \hat {Q}}}_{r \to s} \exp \left( {\xi {\rm {\bf K}}_r } \right),\quad\mbox{(shifting)}$$ in which $\xi$ is the amount of shift over $x$-axis and ${\rm {\bf \hat {Q}}}_{r \to s}$ is the transfer matrix corresponding to shifted space. This property follows (\ref{eq15}) and (\ref{eq21}). The shift theorem plays an important role in the theory of wave propagation in one-dimensional periodic structures in electromagnetic and optics [31]. In the next section we shall show that how one can solve (\ref{eq1}) in its most general form by extension of jump transfer matrices into differential form. \section{Differential Transfer Matrix Method (DTMM)} \subsection{Formulation } First assume that the wavenumber matrix ${\rm {\bf K}}(x)$ is a smoothly varying function of $x$. For the moment, we also assume that $\mathbb{S}$ is non-degenerate. Following the notation in the previous section, let $A$ and $B$ denote respectively the sub-domains $\mathbb{S}_A = \left[ {x - \delta x,x} \right[$ and $\mathbb{S}_B = \left[ {x,x + \delta x} \right[$. If $\delta x$ is small enough, one can approximate the wavenumber matrix as ${\rm {\bf K}}\left( {t \in \mathbb{S}_A } \right) \approx {\rm {\bf K}}(x) \equiv {\rm {\bf K}}_A$, and ${\rm {\bf K}}\left( {t \in \mathbb{S}_B } \right) \approx {\rm {\bf K}}\left( {x + \delta x} \right) \equiv {\rm {\bf K}}_B$. Accordingly, one has $$\label{eq23} {\rm {\bf F}}\left( {x + \delta x} \right) \approx {\rm {\bf F}}_B = {\rm {\bf Q}}_{A \to B} {\rm {\bf F}}_A \approx {\rm {\bf Q}}_{A \to B} {\rm {\bf F}}(x)$$ The above equation is accurate to the first order. So, one can expand the transfer matrix ${\rm {\bf Q}}_{A \to B}$ through its definition (\ref{eq15}) as \label{eq24} \begin{aligned} {\rm {\bf Q}}_{A \to B} & = \exp \left( { - x{\rm {\bf K}}_B } \right)\,{\rm {\bf D}}_B ^{ - 1}{\rm {\bf D}}_A \exp \left( {x{\rm {\bf K}}_A } \right) \\ & = \exp \left[ { - x\left( {{\rm {\bf K}}_A + \delta {\rm {\bf K}}} \right)} \right]\,( {{\rm {\bf D}}_A + \delta {\rm {\bf D}}} )^{ - 1}{\rm {\bf D}}_A \exp \left( {x{\rm {\bf K}}_A } \right) \\ & \approx \exp \left( { - x{\rm {\bf K}}_A } \right)\left( {{\rm {\bf I}} - x\delta {\rm {\bf K}}} \right)\,\left[ {{\rm {\bf D}}_A ^{ - 1} - {\rm {\bf D}}_A ^{ - 1}(\delta {\rm {\bf D}}){\rm {\bf D}}_A ^{ - 1}} \right]{\rm {\bf D}}_A \exp \left( {x{\rm {\bf K}}_A } \right) \\ & = \exp \left( { - x{\rm {\bf K}}_A } \right)\left( {{\rm {\bf I}} - x\delta {\rm {\bf K}}} \right)\,\left( {{\rm {\bf I}} - {\rm {\bf D}}_A ^{ - 1}\delta {\rm {\bf D}}} \right)\exp \left( {x{\rm {\bf K}}_A }\right) \\ & \approx \exp \left( { - x{\rm {\bf K}}_A } \right)\left( {{\rm {\bf I}} - x\delta {\rm {\bf K}} - {\rm {\bf D}}_A ^{ - 1}\delta {\rm {\bf D}}} \right)\exp \left( {x{\rm {\bf K}}_A } \right), \end{aligned} in which we have used the Taylor expansion of matrix exponential and differential property of inverse matrices [32], and have neglected the 2nd- and higher-order variations. The matrix $\delta {\rm {\bf K}}$ can be approximated as $$\label{eq25} \delta {\rm {\bf K}} = \left[ {\delta k_i \delta _{ij} } \right] \approx \left[ {{k}'_i (x)\delta _{ij} } \right] \delta x = \frac{d}{dx}{\rm {\bf K}}(x)\delta x.$$ >From (\ref{eq13}) one also has \label{eq26} \begin{aligned} \delta {\rm {\bf D}} &= \left[ \begin{matrix} 0 & 0 & \cdots & 0 \\ 1 & 1 & \cdots & 1 \\ \vdots & \vdots & & \vdots \\ {\left( {n - 1} \right)_A k_1 ^{n - 2}} & {\left( {n - 1} \right)_A k_2 ^{n - 2}} & \cdots & {\left( {n - 1} \right)_A k_n ^{n -2}} \end{matrix} \right] \frac{d}{dx}{\rm {\bf K}}(x)\delta x \\ &\equiv {\rm {\bf C}}\delta {\rm {\bf K}}. \end{aligned} So we can rewrite (\ref{eq23}) as $$\label{eq27} \delta {\rm {\bf F}}(x) \approx \left( {{\rm {\bf Q}}_{A \to B} - {\rm {\bf I}}} \right){\rm {\bf F}}(x) \equiv {\rm {\bf U}}(x){\rm {\bf F}}(x)\delta x.$$ But the matrix ${\rm {\bf U}}(x)$ can be found from (\ref{eq24}) and (\ref{eq26}) as $$\label{eq28} {\rm {\bf U}}(x) \approx - \exp \left( { - x{\rm {\bf K}}} \right)\left( {x{\rm {\bf I}} + {\rm {\bf D}}^{ - 1}{\rm {\bf C}}} \right)\frac{\delta {\rm {\bf K}}}{\delta x }\exp \left( {x{\rm {\bf K}}} \right) ,$$ where the trivial subscript $A$ has been dropped. Now, if we let $\delta x$ approach zero, (\ref{eq27}) and (\ref{eq28}) become exact and can be rewritten as $$\label{eq29} d{\rm {\bf F}}(x) = {\rm {\bf U}}(x){\rm {\bf F}}(x)dx,$$ $$\label{eq30} {\rm {\bf U}}(x) = - x{\rm {\bf {K}'}}(x) - \exp \left[ { - x{\rm {\bf K}}(x)} \right]{\rm {\bf D}}(x)^{ - 1}{\rm {\bf C}}(x){\rm {\bf {K}'}}(x)\exp \left[ {x{\rm {\bf K}}(x)} \right].$$ Hereinafter, we refer to ${\rm {\bf U}}(x)$ as the differential transfer or the kernel matrix. In (\ref{eq30}), we have \begin{gather} \label{eq32} {\rm {\bf C}}(x) = \left[ {\left( {i - 1} \right)k_j^{i - 2} (x)} \right]_{n\times n}, {\rm {\bf D}}(x) = \left[ {k_j^{i - 1} (x)} \right]_{n\times n}, \\ {\rm {\bf K}}(x) = \left[ {k_i (x)\delta _{ij} } \right]_{n\times n}, {\rm {\bf {K}'}}(x) = \left[ {{k}'_i (x)\delta _{ij} } \right]_{n\times n}. \end{gather} While (\ref{eq29}) can be integrated by numerical methods [23-25], it permits an exact solution of the form $$\label{eq310} {\rm {\bf F}}(x_2) = {\rm {\bf Q}}_{x_1 \to x_2 } {\rm {\bf F}}(x_1),$$ for $\forall x_1 ,x_2 \in \mathbb{S}.$. Here ${\rm {\bf Q}}_{x_1 \to x_2 }$ is referred to as the transfer matrix from $x_1$ to $x_2$, given by the Dyson's perturbation theory as [33-37] \label{eq311} \begin{aligned} {\rm {\bf Q}}_{x_1 \to x_2 } &={\bf I} + \int_{x_1 }^{x_2}dt_1{\rm {\bf U}}(t_1) + \int_{x_1 }^{x_2}dt_1 \int_{x_1 }^{t_1}dt_2{\rm {\bf U}}(t_1){\rm {\bf U}}(t_2) + \ldots \\ &\quad + \int_{x_1 }^{x_2}dt_1 \int_{x_1 }^{t_1}dt_2 \int_{x_1 }^{t_2}dt_3 \ldots \int_{x_1 }^{t_{m-1}}dt_m {\rm {\bf U}}(t_1) \ldots {\rm {\bf U}}(t_m) + \ldots, \end{aligned} By defining $\mathbb{P}$ as the so-called chronological ordering operator [33-37], it is possible to rewrite (\ref{eq311}) as $$\label{eq311a} {\rm {\bf Q}}_{x_1 \to x_2 } = \sum_{m=0}^{\infty}\frac{1}{m!} \int_{x_1 }^{x_2}dt_1 \int_{x_1 }^{x_2}dt_2 \int_{x_1 }^{x_2}dt_3 \ldots \int_{x_1 }^{x_2}dt_m \mathbb{P} \left[ {\rm {\bf U}}(t_1) \ldots {\rm {\bf U}}(t_m) \right].$$ Often, (\ref{eq311a}) is symbolically written as [33-37] $$\label{eq31} {\rm {\bf Q}}_{x_1 \to x_2 } = \mathbb{P} \exp \Big[ {\int_{x_1 }^{x_2 } {{\rm {\bf U}}(x)dx} } \Big] \equiv \mathbb{P} \exp \left( {{\rm {\bf M}}_{x_1 \to x_2 } } \right),$$ in which ${\rm {\bf M}}_{x_1 \to x_2}$ is the integral of the kernel matrix and referred to as the transfer exponent matrix. The above expression can be greatly simplified if the kernel matrix ${\rm {\bf U}}(x)$ satisfies one of the few existing sufficient conditions for integrability, including when it commutes with ${\rm {\bf M}}_{x_1 \to x_2}$ [38-40], or it satisfies the Fedorov's condition [40], or it commutes with itself for $\forall x_1, x_2 \in \mathbb{S}$. In this case, one has $$\label{eq313} {\rm {\bf Q}}_{x_1 \to x_2 } = \exp \Big[ {\int_{x_1 }^{x_2 } {{\rm {\bf U}}(x)dx} } \Big].$$ Although in general (\ref{eq313}) is not necessarily the exact transfer matrix, it is known that both (\ref{eq310}) and (\ref{eq313}) must have the same determinant [39]. Also their traces are equal at least to the second order. This can be easily justified by comparing the expansion of (\ref{eq313}) to (\ref{eq311}) and using the cyclic permutation property of $\mathop{\rm tr}\left\{\cdot\right\}$ [30]. This means for the case of second-order equations with $n=2$, (\ref{eq313}) can be used instead of (\ref{eq31}), with an accuracy better than second-order. If needed, the exact transfer matrix can be directly found by numerical integration of equation [41] $$\label{eq314} d{\rm {\bf Q}}_{x_1 \to x } = {\rm {\bf U}}(x) {\rm {\bf Q}}_{x_1 \to x } dx,$$ with the initial condition ${\rm {\bf Q}}_{x_1 \to x_1}={\bf I}$. \subsection{Properties of Transfer Matrix} The transfer matrix ${\rm {\bf Q}}_{x_1 \to x_2 }$ from $x_1$ to $x_2$ clearly preserves the properties (\ref{eq21a}), (\ref{eq21b}) and (\ref{eq21c}), and as discussed above, its determinant is always equal to the determinant of its counterpart given by (\ref{eq313}). Therefore, $\left|{\rm {\bf Q}}_{x_1 \to x_2}\right|=\exp(\mathop{\rm tr}\left\{{\rm {\bf M}}_{x_1 \to x_2 } \right\})$. But the transfer exponent matrix ${\rm {\bf M}}_{x_1 \to x_2 }$ can be written as the sum of the jump ${\rm {\bf J}}_{x_1 \to x_2 }$ and propagation ${\rm {\bf T}}_{x_1 \to x_2 }$ matrices given by \begin{gather} \label{eq33} {\rm {\bf J}}_{x_1 \to x_2 } = - \int_{x_1 }^{x_2 } {\exp \left[ { - x{\rm {\bf K}}(x)} \right]{\rm {\bf D}}(x)^{ - 1}\frac{d}{dx}\left\{ {{\rm {\bf D}}(x)\exp \left[ {x{\rm {\bf K}}(x)} \right]} \right\}\,dx},\\ \label{eq34} {\rm {\bf T}}_{x_1 \to x_2 } = \int_{x_1 }^{x_2 } {{\rm {\bf K}}(x)dx}. \end{gather} The validity of the identity ${\rm {\bf M}}_{x_1 \to x_2 }={\rm {\bf J}}_{x_1 \to x_2 }+{\rm {\bf T}}_{x_1 \to x_2 }$ can be verified by comparing the integrands of both sides. The propagation matrix ${\rm {\bf T}}_{x_1 \to x_2 }$ is diagonal and also has a diagonal matrix exponential given by $$\label{eq35} \exp \left( {{\rm {\bf T}}_{x_1 \to x_2 } } \right) = \Big[ {\exp \Big( {\int_{x_1 }^{x_2 } {k_i (x)dx} } \Big)\delta _{ij} } \Big].$$ Therefore, its determinant can be easily found by multiplying its diagonal elements. The determinant $\left|{\rm {\bf J}}_{x_1 \to x_2 }\right|$ can be obtained using the relation ${\rm {\bf {D}'}}(x) = {\rm {\bf C}}(x){\rm {\bf {K}'}}(x)$ according to (\ref{eq26}), and the identity (see Appendix A) $$\label{eq36} \Big| \exp \Big[ {\int_{x_1}^{x_2} {{\rm {\bf H}}(x)^{ - 1}{\rm {\bf {H}'}}(x)dx} } \Big] \Big| = \big|{\rm {\bf H}}^{ - 1}\left( x_1 \right){\rm {\bf H}}\left( x_2 \right) \big|,$$ as $$\label{eq37} \big|\exp \left( {{\rm {\bf J}}_{x_1 \to x_2 } } \right) \big|= \big| \exp \left[ { - x_2 {\rm {\bf K}}(x_2)} \right]{\rm {\bf D}}(x_2)^{ - 1}{\rm {\bf D}}(x_1)\exp \left[ {x_1 {\rm {\bf K}}(x_1)} \right] \big|.$$ Finally, the determinant $\left| {\rm {\bf Q}}_{x_1 \to x_2 }\right|$ can be thus found using the well-known identities [29] $\left| {\exp \left( {{\rm {\bf A}} + {\rm {\bf B}}} \right)} \right| = \left| {\exp \left( {\rm {\bf A}} \right)} \right|\times \left| {\exp \left( {\rm {\bf B}} \right)} \right|$, and $\left| {\exp \left( {\rm {\bf M}} \right)} \right| = \exp \left( {\mathop{\rm tr} \left\{ {\rm {\bf M}} \right\}} \right)$ as $$\label{eq38} \left| {{\rm {\bf Q}}_{x_1 \to x_2 } } \right| = \exp \left( {\mathop{\rm tr}\left\{ {{\rm {\bf T}}_{x_1 \to x_2 } } \right\}} \right)\times \exp \left( {\mathop{\rm tr}\left\{ {{\rm {\bf J}}_{x_1 \to x_2 } } \right\}} \right).$$ Using the (\ref{eq35}), (\ref{eq37}), and (\ref{eq17}) one can finally obtain \label{eq39} \begin{aligned} &\left| {{\rm {\bf Q}}_{x_1 \to x_2 } } \right| \\ &= \exp \Big( {x_1 \mathop{\rm tr}\left\{ {{\rm {\bf K}}(x_1)} \right\} - x_2 \mathop{\rm tr}\left\{ {{\rm {\bf K}}(x_2)} \right\} - \int_{x_1 }^{x_2 } {\mathop{\rm tr}\left\{ {{\rm {\bf K}}(x)} \right\}dx} } \Big)\times \prod_{i > j} {\frac{k_i (x_1) - k_j (x_1)}{k_i (x_2) - k_j (x_2)}} \\ &= \exp \Big[ {x_1 a_{n - 1} (x_1) - x_2 a_{n - 1} (x_2) - \int_{x_1 }^{x_2 } {a_{n - 1} (x)dx} } \Big]\times \prod_{i > j} {\frac{k_i (x_1) - k_j (x_1)}{k_i (x_2) -k_j (x_2)}} .\\ \end{aligned} This equation can be compared to (\ref{eq21d}) as the determinant property of transfer matrix. Notice that if $a_{n - 1} (x) \equiv 0$, then (\ref{eq39}) reduces to $$\label{eq40} \left| {{\rm {\bf Q}}_{x_1 \to x_2 } } \right| = \prod_{i > j} {\frac{k_i (x_1) - k_j (x_1)}{k_i (x_2) - k_j (x_2)}}.$$ This shows that the determinant of the transfer matrix is only a function of starting and end points $x_1$ and $x_2$. Note that the transformation $$\label{eq41} f(x) \to \exp \Big[ { - {1 \over n}\int {a_{n - 1} (x)dx} } \Big] h(x)$$ in (\ref{eq1}) always gives rise to a new $n$th-order ODE with identically vanishing coefficient of $h^{\left( {n - 1} \right)}(x)$ term. Therefore, the above property of determinants can be always met through the above mentioned transformation. The shifting property (\ref{eq22a}) and scaling properties must also preserve, since the transfer matrix ${\rm {\bf Q}}_{x_1 \to x_2 }$is indeed obtained by dividing the sub-domain $\left[ {x_1 ,x_2 } \right]$ into infinitely many thin layers and multiplying the infinitesimal jump transfer matrices. \subsection{Justification of DTMM } Here we prove that the DTMM formulation is exact by showing that the solution satisfies (\ref{eq1}) through direct substitution. \begin{lemma} \label{lm1} Let the Vandermonde matrix \textbf{D} be invertible. Then the elements of its inverse ${\rm {\bf D}}^{ - 1} = \left[ {\gamma _{ij} } \right]$ satisfy $\sum_{i = 1}^n {\gamma _{ir} k_i ^{m - 1}} = \delta _{mr}$ where $1 \le m \le n$. \end{lemma} \begin{proof} The Vandermonde coefficients $\gamma _{ij}$ are found from the expansion of Lagrange interpolation polynomials as [42,43] $$\label{eq42} \Gamma _i (t) = \prod_{j \ne i} {\frac{t - k_j }{k_i - k_j }} = \sum_{j = 1}^n {\gamma _{ij} t^{j - 1}} ,\quad i = 1,\ldots, n.$$ Obviously $\Gamma _i \left( {k_j } \right) = \sum_{r = 1}^n {\gamma _{ir} k_j ^{r - 1}} = \delta _{ij}$. Therefore $k_j ^{m - 1} = \sum_{i = 1}^n {k_i ^{m - 1}\Gamma _i \left( {k_j } \right)}$, or $$\label{eq43} k_j ^{m - 1} = \sum_{i = 1}^n {k_i ^{m - 1}\sum_{r = 1}^n {\gamma _{ir} k_j ^{r - 1}} } = \sum_{r = 1}^n {\Big( {\sum_{i = 1}^n {\gamma _{ir} k_i ^{m - 1}} } \Big)k_j ^{r - 1}} \equiv \sum_{r = 1}^n {a_{mr} k_j ^{r - 1}}.$$ Given a fixed $m$, the above equation can be transformed into a linear system of equations as ${\rm {\bf A}}^t{\rm {\bf D}} = {\rm {\bf B}}^t$, with \textbf{D} being the Vandermonde matrix, ${\rm {\bf A}} = \left[ {a_{rm} } \right]$ being the column vector of unkowns, and ${\rm {\bf B}} = \left[ {k_i ^{m - 1}} \right]$ being the input column vector. Since \textbf{D} is invertible by assumption, then \textbf{A} and thus $a_{rm}$ must be unique. However, (\ref{eq43}) is satisfied through setting $a_{rm} = \delta _{rm}$ for $1 \le m \le n$, and hence the proof is complete. \end{proof} \begin{lemma} \label{lm2} Suppose that ${\rm{\bf F}}(x)$ is an $n\times 1$ vector function, defined on the connected non-degenerate domain $\mathbb{S}$, satisfying the differential equation \eqref{eq29} with the kernel matrix given by \eqref{eq30}. Then, we have $$\label{eq44} \frac{d^m}{dx^m}\left\{ {\exp \left[ {{\rm {\bf \Phi }}(x)} \right]^{\,t}{\rm {\bf F}}(x)} \right\} = \exp \left[ {{\rm {\bf \Phi }}(x)} \right]^{\,t}{\rm {\bf K}}^m(x){\rm {\bf F}}(x), m = 0,\ldots, n,$$ where ${\rm {\bf K}}(x)$ is the wavenumber matrix given by \eqref{eq9} and ${\rm {\bf\Phi }}(x)$ is the corresponding phase vector defined in \eqref{eq7}. \end{lemma} \begin{proof} For $m = 0$ the proof is trivial. Also for $m \ge 0$ we have \label{eq45} \begin{aligned} & \frac{d}{dx}\left\{ {\exp \left( {\rm {\bf \Phi }} \right)^{\,t}{\rm {\bf K}}^m{\rm {\bf F}}} \right\} \\ &= \frac{d}{dx}\exp \left( {\rm {\bf \Phi }} \right)^{\,t}{\rm {\bf K}}^m{\rm {\bf F}} + \exp \left( {\rm {\bf \Phi }} \right)^{\,t}\frac{d}{dx}{\rm {\bf K}}^m{\rm {\bf F}} + \exp \left( {\rm {\bf \Phi }} \right)^{\,t}{\rm {\bf K}}^m\frac{d}{dx}{\rm {\bf F}} \\ & = \exp \left( {\rm {\bf \Phi }} \right)^{\,t}\left( {{\rm {\bf K}}^{m + 1} + x{\rm {\bf K}}^m{\rm {\bf {K}'}} + m{\rm {\bf K}}^{m - 1}{\rm {\bf {K}'}} + {\rm {\bf K}}^m{\rm {\bf U}}} \right){\rm {\bf F}}\,, \end{aligned} where the dependence of matrices on $x$ is not shown for the sake of convenience. From the definition of the kernel matrix (\ref{eq30}), (\ref{eq45}) can be simplified as \label{eq46} \begin{aligned} &\frac{d}{dx}\left\{ {\,\exp \left( {\rm {\bf \Phi }} \right)^{\,t}{\rm {\bf K}}^m{\rm {\bf F}}} \right\} \\ &= \exp \left( {\rm {\bf \Phi }} \right)^{\,t}{\rm {\bf K}}^{m+1}{\rm {\bf F}} \\ &\quad + \exp \left( {\rm {\bf \Phi }} \right)^{\,t}\left[ {m{\rm {\bf K}}^{m - 1}{\rm {\bf {K}'}} - {\rm {\bf K}}^m\exp \left( { - x{\rm {\bf K}}} \right){\rm {\bf D}}^{ - 1}{\rm {\bf C{K}'}}\exp \left( {x{\rm {\bf K}}} \right)} \right]\,{\rm {\bf F}} \\ & = \exp \left( {\rm {\bf \Phi }} \right)^{\,t}{\rm {\bf K}}^{m + 1}{\rm {\bf F}} + {\rm {\bf 1}}_{n\times 1}^{\,t} \left( {m{\rm {\bf K}}^{m - 1} - {\rm {\bf K}}^m{\rm {\bf D}}^{ - 1}{\rm {\bf C}}} \right)\,{\rm {\bf {K}'}}\exp \left( {x{\rm {\bf K}}} \right){\rm {\bf F}} \\ & \equiv \exp \left( {\rm {\bf \Phi }} \right)^{\,t}{\rm {\bf K}}^{m + 1}{\rm {\bf F}} + {\rm {\bf 1}}_{1\times n} {\rm {\bf V{K}'}}\exp \left( {x{\rm {\bf K}}} \right){\rm {\bf F}} \,, \end{aligned} where ${\rm {\bf 1}}_{n\times 1}$ is defined under (\ref{eq9}) and ${\rm {\bf V} }\equiv {m{\rm {\bf K}}^{m - 1} - {\rm {\bf K}}^m{\rm {\bf D}}^{ - 1}{\rm {\bf C}}}$ . Now we show that the row vector ${\rm {\bf 1}}_{1\times n} {\rm {\bf V}}$ is identically zero, but only if $m < n$. By defining the elements of the inverse of the Vandermonde matrix as ${\rm {\bf D}}^{ - 1} = \left[ {\gamma _{ij} } \right]$, we have \label{eq47} \begin{aligned} {\rm {\bf 1}}_{1\times n} {\rm {\bf V}} &= m{\rm {\bf 1}}_{1\times n} {\rm {\bf K}}^{m - 1} - {\rm {\bf 1}}_{1\times n} {\rm {\bf K}}^m{\rm {\bf D}}^{ - 1}{\rm {\bf C}} \\ & = \Big[ {\sum_{i =1}^n {mk_i ^{m - 1}\delta _{ij} } } \Big] - \Big[ {\sum_{i =1}^n {\sum_{r = 1}^n {\gamma _{ir} \left( {r - 1} \right)k_j ^{r - 2}k_i ^m} } } \Big] \\ & = \Big[ {mk_j ^{m - 1} - \sum_{r = 1}^{n - 1} {rk_j ^{r - 1}\sum_{i = 1}^n {\gamma _{ir + 1} k_i ^m} } } \Big] \,. \end{aligned} Since the inverse of \textbf{D} must exist by assumption, we thus have according to Lemma (\ref{lm1}) $$\label{eq48} {\rm {\bf 1}}_{1\times n} {\rm {\bf V}} = \Big[ {mk_j ^{m - 1} - \sum_{r = 1}^{n - 1} {rk_j ^{r - 1}\delta _{mr} } } \Big] = {\rm {\bf 0}}_{1\times n},\; m < n .$$ Therefore, the proof of theorem follows immediately by induction. \end{proof} \begin{theorem} [Fundamental theorem of differential transfer matrix method] \label{thm1} Suppose that ${\rm {\bf F}}(x)$ satisfies the derivative Lemma, with the kernel matrix ${\rm {\bf K}}(x)$ satisfying the characteristic equation \eqref{eq8}. Then, the function $f(x)$ defined by \eqref{eq6} is a solution of the differential equation \eqref{eq1}. \end{theorem} \begin{proof} By the expansion of the operator $\mathbb{L}$ given in (\ref{eq2}) on $f(x)$ given in (\ref{eq6}), and using the derivative Lemma we have \label{eq49} \begin{aligned} {\rm \mathbb{L}}f(x) &= \sum_{m = 0}^n {a_m (x)\frac{d^m}{dx^m}} f(x)\\ &= \sum_{m = 0}^n {a_m (x)\frac{d^m}{dx^m}\left\{ {\exp \left[ {{\rm {\bf \Phi }}(x)} \right]^{\,t}{\rm {\bf F}}(x)} \right\}} \\ &= \sum_{m = 0}^n {a_m (x)\exp \left[ {{\rm {\bf \Phi }}(x)} \right]^{\,t}{\rm {\bf K}}^m(x){\rm {\bf F}}(x)}\\ &= \exp \left[ {{\rm {\bf \Phi }}(x)} \right]^{\,t}\Big[ {\sum_{m = 0}^n {a_m (x){\rm {\bf K}}^m(x)} } \Big]{\rm {\bf F}}(x) \,. \end{aligned} However, the summation within the brackets is the null matrix because of (\ref{eq8}). Thus the right-hand-side of the above equation must be identically zero. This completes the proof. \end{proof} \subsection{Linear Independent Solutions} Consider the set of vectors ${\rm {\bf F}}_i \in \mathbb{S}^n,i = 1,\ldots, n$, forming a basis in $\mathbb{S}^n$. If ${\bf Q}$ is a non-singular matrix, then the vectors ${\rm {\bf G}}_i = {\rm {\bf QF}}_i \in \mathbb{C}^n,i = 1,\ldots, n$ would clearly constitute another basis on $\mathbb{S}^n$. Let $\mathbb{S}$ be a non-degenerate domain with the differential equation (\ref{eq1}) and solution given by (\ref{eq310}). Since $\mathbb{S}$ is non-degenerate, then ${\rm {\bf Q}}_{x_1 \to x_2 }$ must be non-singular for all $x_1 ,x_2 \in \mathbb{S}$. Therefore the set of $n$ vector functions defined by ${\rm {\bf G}}_i (x) = {\rm {\bf Q}}_{a \to x} {\rm {\bf F}}_i \left( a \right) \in \mathbb{S}^n$ form a basis on $\mathbb{S}^n$. Let $\Psi$ denote the set of all functions operating on $\mathbb{S}$, as $\Psi = \left\{ {\psi :\mathbb{S} \mapsto \mathbb{C}} \right\}$. We define the scalar functions $g_i \in \Psi$ as $g_i (x) = \exp \left[ {{\rm {\bf \Phi }}(x)} \right]^{\,t}{\rm {\bf G}}_i (x)$ and show that they form a basis, by inspection of their Wronskian determinant. \begin{theorem} \label{thm2} With the above definitions, the set of functions $g_i (x), i = 1,\ldots, n$ form $n$ linear independent solutions on $\Psi$. \end{theorem} \begin{proof} The Wronskian matrix in $\Psi$ can be written as ${\rm {\bf W}} = \left[ {g_i ^{(j-1)}(x)} \right]$. Using Theorem 1 one has ${\rm {\bf W}} = \left[ {\exp \left( {\rm {\bf \Phi }} \right){\rm {\bf K}}^{j - 1}{\rm {\bf Q}}_{a \to x} {\rm {\bf F}}_i } \right]$, with the dependences on $x$ omitted. Now we define the $n\times n$ matrix ${\rm {\bf F}} = \left[ {{\rm {\bf F}}_i } \right]$. Then ${\rm {\bf W}} = {\rm {\bf D}}\exp \left( {x{\rm {\bf K}}} \right)\,{\rm {\bf Q}}_{a \to x} \,{\rm {\bf F}}$. Because $\mathbb{S}$ is non-degenerate $\left| {\rm {\bf D}} \right| \ne 0$. Also, ${\rm {\bf Q}}_{a \to x}$ must be non-singular and thus $\left| {\rm {\bf Q}}_{a \to x}\right| \ne 0$. Now since according to the assumption, ${\rm {\bf F}}_i$ form a basis on $\mathbb{S}^n$, $\left| {\rm {\bf F}} \right| \ne 0$ and therefore for the Wronskian determinant we have $\left| {\rm {\bf W}} \right| \ne 0$. Thus, the proof is complete. \end{proof} \subsection{Treatment of Singularities} Now let $\mathbb{S}$ is degenerate at a finite number of isolated singularities $\xi _i \in \mathbb{S},i = 1,\ldots, \Sigma$, at which $\left| {{\rm {\bf D}}\left( {\xi _i } \right)} \right| = 0$. Therefore, the ${\rm {\bf U}}\left( {\xi _i } \right)$ is singular and ${\rm {\bf M}}_{x_1 \to x_2 }$ does not exist. Without loss of generality, we can assume $\Sigma = 1$. If $\mathbb{S}$ represents the integration domain $\left[ {x_1 ,x_2 } \right]$, the total transfer matrix ${\rm {\bf Q}}_{x_1 \to x_2 }$ by the decomposition property (\ref{eq21c}) can be written as $$\label{eq50} {\rm {\bf Q}}_{x_1 \to x_2 } = {\rm {\bf Q}}_{\xi + \delta x \to x_2 } {\rm {\bf Q}}_{\xi - \delta x \to \xi + \delta x} {\rm {\bf Q}}_{x_1 \to \xi - \delta x}.$$ The transfer matrices ${\rm {\bf Q}}_{x_1 \to \xi - \delta x}$ and ${\rm {\bf Q}}_{\xi + \delta x \to x_2 }$ do exist and can be evaluated directly by integration and exponentiation of ${\rm {\bf U}}(x)$, as long as $\delta x$ is positive and finite. The transfer matrix ${\rm {\bf Q}}_{\xi - \delta x \to \xi + \delta x}$ enclosing the singularity, also can be approximated by its equivalent jump transfer matrix as \label{eq51} \begin{aligned} &{\rm {\bf Q}}_{\xi - \delta x \to \xi + \delta x}\\ & \approx \exp \left[ { - \left( {\xi + \delta x} \right){\rm {\bf K}}\left( {\xi + \delta x} \right)} \right]\,{\rm {\bf D}}^{ - 1}\left( {\xi + \delta x} \right) {\rm {\bf D}}\left( {\xi - \delta x} \right)\exp \left[ {\left( {\xi - \delta x} \right){\rm {\bf K}}\left( {\xi- \delta x} \right)} \right]\,. \end{aligned} This approach permits evaluation of the total transfer matrix ${\rm {\bf Q}}_{x_1 \to x_2 }$ by making finite jumps over singularities. If $\mathbb{S}$ is entirely degenerate, then the transformation $f(x) \to w(x)h(x)$ with $w(x)$ being a non-constant function, results in a completely different characteristic equation (\ref{eq8}). A proper choice of $w(x)$ leads to a new $\mathbb{S}$, which can be no longer entirely degenerate. Then DTMM can be used to solve for $h(x)$. \section{Examples} \subsection{First-order ODEs} It is easy to show the consistency of the method for first order differential equations with variable coefficients, that is $$\label{eq55} {f}'(x) + a_0 (x)f(x) = 0,$$ having the exact solution $$\label{eq56} f(x_2) = \exp \Big( {\int_{x_2 }^{x_1 } {a_0 (x)dx} } \Big)f(x_1).$$ In this case, the only root of (\ref{eq8}) is $k(x)=-a_0(x)$, all matrices become scalar with $D = 1$ and $C = 0$, and the transfer exponent reduces to $$\label{eq57} M_{x_1 \to x_2 } = \int_{x_1 }^{x_2 } {x\frac{d}{dx}\left[ {a_0 (x)} \right]dx} = x_2 a_0 (x_2) - x_1 a_0 (x_1) - \int_{x_1 }^{x_2 } {a_0 (x)dx}.$$ Clearly (\ref{eq313}) holds and therefore, accordingly one has $$\label{eq58} F(x_2) = \exp \Big[ {x_2 a_0 (x_2) - x_1 a_0 (x_1) - \int_{x_1 }^{x_2 } {a_0 (x)dx} } \Big]F(x_1).$$ From (\ref{eq6}), we have $$\label{eq59} f(x) = \exp \left[ {xk(x)} \right] F(x).$$ Therefore, using (\ref{eq58}) and (\ref{eq59}) one can obtain (\ref{eq56}). \subsection{Second-order ODEs} In this section, we first consider the simplest second-order differential equation with variable coefficients, given by $$\label{eq60} {f}''(x) + a_0 (x)f(x) = 0.$$ This equation has been studied since 1836 by Sturm [44] and has been considered in many reports [45-50]. Actually, any second-order ODE can be rewritten in the form of (\ref{eq60}) by the transformation given in (\ref{eq41}). Also, the first order non-linear Riccati equation [2] takes the above form after suitable transformations. In physics this equation models the propagation of transverse-electric (TE) polarized light in one-dimensional (1D) non-homogeneous media (with a position-dependent refractive index), or motion of a single electron in a 1D potential trap. In optical problems $f(x)$ is the amplitude of transverse electric field and $a_0 (x) = k_0 ^2\epsilon _r (x) - \beta ^2$, in which $k_0$ is the wavenumber of the radiation in vacuum, $\epsilon _r (x)$ is the relative permittivity function of the medium, and $\beta$ is the propagation eigenvalue. In quantum mechanics $f(x)$ is the probability wave function and $a_0 (x) = {2m\left[ {E - V(x)} \right]} \mathord{\left/ {\vphantom {{2m\left[ {E - V(x)} \right]} {\hbar ^2}}} \right. \kern-\nulldelimiterspace} {\hbar ^2}$, where $m$ is the electron mass, $V(x)$ is the electric potential, $\hbar$ is the reduced Planck constant, and $E$ is the energy. Radial functions of axial gravity waves of a black hole are governed by a generalized Klein-Gordon or Regge-Wheeler equation [51,52], which takes a similar form to (\ref{eq60}). Here, we consider the solution of this equation by DTMM. For $n = 2$, the kernel matrix ${\rm {\bf U}}(x)$ can be simplified from (\ref{eq30}) to $$\label{eq61} {\rm {\bf U}}(x) = \left[ \begin{matrix} { - \left( {x + \frac{1}{k_1 - k_2 }} \right){k}'_1 } & {\frac{{k}'_2 }{k_2 - k_1 }\exp \left[ { - x\left( {k_1 - k_2 } \right)} \right]} \\ {\frac{{k}'_1 }{k_1 - k_2 }\exp \left[ { + x\left( {k_1 - k_2 } \right)} \right]} & { - \big( {x + \frac{1}{k_2 - k_1 }} \big){k}'_2 } \end{matrix}\right].$$ From (\ref{eq8}), one has $k_i^2(x)+a_0(x)=0,i=1,2$, that is $k_1 (x) = - jk(x)$ and $k_2 (x) = + jk(x)$, with $k(x) \equiv \sqrt {a_0 (x)}$. Therefore (\ref{eq61}) simplifies as $$\label{eq62} {\rm {\bf U}}(x) = \frac{{k}'(x)}{2k(x)}\left[ \begin{matrix} { - 1 + j2k(x)x} & {\exp \left[ {j2xk(x)} \right]} \\ {\exp \left[ { - 2jxk(x)} \right]} & { - 1 - j2k(x)x} \end{matrix}\right],$$ which is in agreement with our previous results [18,19,22]. The transfer matrix ${{\rm {\bf Q}}_{x_1 \to x_2 } }$ is then found by (\ref{eq31}) (exact evaluation of matrix exponential in this case is possible by the theorem described in the Appendix B). The determinant of ${{\rm {\bf Q}}_{x_1 \to x_2} }$ is $$\label{eq63} \left| {{\rm {\bf Q}}_{x_1 \to x_2 } } \right| = \frac{k_2 (x_1) - k_1 (x_1)}{k_2 (x_2) - k_1 (x_2)} = \frac{k(x_1)}{k(x_2)}.$$ The complete solution of (\ref{eq60}) is then given by (\ref{eq5}) as $$\label{eq64} f(x) = f_1 (x)\exp \left[ { - jxk(x)} \right] + f_2 (x)\exp \left[ { + jxk(x)} \right].$$ The singularities of (\ref{eq60}) correspond to the turning points at which both wavenumbers become zero. In fact, a singularity separates the regions in which the waves are evanescent and propagating, or in other words, $a_0 (x)$ is respectively positive or negative. The treatment method for singularities as described in sub-section 4.5 can be applied to find the transfer matrix over the singularity, being approximated by its corresponding jump transfer matrix. We can classify the singularities depending on the nature of the waves across the singularity. Here, a singularity at $x = \xi$ is characterized as $k_1(\xi)=k_2(\xi)$, or $a_0(\xi)=0$. Correspondingly, the singularity at $x = \xi$ is referred to as type A if $a_0 \left( {x > \xi } \right) > 0$ and $a_0 \left( {x < \xi } \right) < 0$, type B if $a_0 \left( {x > \xi } \right) < 0$ and $a_0 \left( {x < \xi } \right) > 0$, and otherwise type C. Due to (\ref{eq15}) the jump transfer matrix from region 1 to 2 across the interface $x = \xi$ for (\ref{eq60}) takes the form [18,19,22,53] $$\label{eq65} {\rm {\bf Q}}_{1 \to 2} = \left[ \begin{matrix} {\frac{k_2 + k_1 }{2k_2 }e^{ + j\xi \left( {k_2 - k_1 } \right)}} & {\frac{k_2 - k_1 }{2k_2 }e^{ + j\xi \left( {k_2 + k_1 } \right)}} \\ {\frac{k_2 - k_1 }{2k_2 }e^{ - j\xi \left( {k_2 + k_1 } \right)}} & {\frac{k_2 + k_1 }{2k_2 }e^{ - j\xi \left( {k_2 - k_1 } \right)}} \\ \end{matrix}\right].$$ Here, $k_1^2=-a_0(\xi-\delta x)$ and $k_1^2=-a_0(\xi+\delta x)$. Therefore, the jump transfer matrix over singularities given by ${\rm {\bf Q}}_{\xi - \delta x \to \xi + \delta x}$ simplifies as \begin{gather} \label{eq66a} {\rm {\bf Q}}_{\xi - \delta x \to \xi + \delta x} = \left[ \begin{matrix} {\frac{1 + j}{2}} & {\frac{1 - j}{2}} \\ {\frac{1 - j}{2}} & {\frac{1 + j}{2}} \end{matrix} \right],\quad\mbox{type A}\\ \label{eq66b} {\rm {\bf Q}}_{\xi - \delta x \to \xi + \delta x} = \left[ \begin{matrix} {\frac{1 - j}{2}} & {\frac{1 + j}{2}} \\ {\frac{1 + j}{2}} & {\frac{1 - j}{2}} \end{matrix}\right], \quad \mbox{type B}\\ \label{eq66c}{\rm {\bf Q}}_{\xi - \delta x \to \xi + \delta x} = \left[ \begin{matrix} 1 & 0 \\ 0 & 1 \\ \end{matrix} \right],\quad \mbox{type C}. \end{gather} The results here are obtained by setting $k\left( {\xi \pm \delta x} \right) \approx \sqrt {\pm \delta x{a}'_0 \left( {\xi \pm \delta x} \right)}$ through the corresponding Taylor's expansion of $a_0 (x)$ around $x = \xi$ and taking the limit $\delta x \to 0^ +$. The DTMM has been used to solve (\ref{eq60}) numerically [18,19], and the results have been compared to other approaches. In cases that analytical solutions are known, the results have been in agreement. While the exact transfer matrix is given by (\ref{eq31}), we have noticed that the reduced form as in (\ref{eq313}) can be also used and leads to accurate solutions. Also, energy dispersion of an electron in the 1D potential of infinite monatomic chain has been computed [22] by this approach. In general, computation times for DTMM are at least one order of magnitude lower than the other existing approaches. We also have discussed how it is possible to extend this approach to transverse-magnetic (TM) polarized modes and inhomogeneous anisotropic media, where four field components are solved for at once [22]. We also have recently exploited the DTMM to periodic electromagnetic structures [39], in which $a_0(x)$, and thus the wavenumber functions $k_i(x),i=1,2$ are all periodic. This has led us to a novel yet simple mathematical description of waves in inhomogeneous periodic media. \subsection{Fourth-order ODEs} As another application example, we obtain the differential transfer matrix of a simple fourth-order differential equation that has been discussed in several reports [49,50] $$\label{eq67} f^{\left( 4 \right)}(x) + a_0 (x)f(x) = 0.$$ This time from (\ref{eq8}) we have $k_i^4(x)+a_0(x)=0, i=1,\ldots,4$, and thus $k_1 (x) = - k(x)$, $k_2 (x) = - jk(x)$, $k_3 (x) = + k(x)$, and $k_4 (x) = + jk(x)$, in which $k(x) \equiv \sqrt[4]{ - a_0 (x)}$. The kernel matrix ${\rm {\bf U}}(x)$ is found by {\tt Mathematica} as \label{eq68} \begin{aligned} &{\rm {\bf U}}(x) =\frac{{k}'(x)}{2k(x)} \times \\ &\left[ \begin{matrix} {2xk(x) - 3} & {(1+j)e^{(1-j)xk(x)}} & {e^{2xk(x)}} & {(1-j)e^{(1+j)xk(x)}} \\ {(1-j)e^{(j-1)xk(x)}} & {2jxk(x) - 3} & {(1+j)e^{(1+j)xk(x)}} & {e^{2jxk(x)}} \\ {e^{ - 2xk(x)}} & {(1-j)e^{ - (1+j)xk(x)}} & { - 2xk(x) - 3} & {(1+j)e^{(j-1)xk(x)}} \\ {(1+j)e^{ - (1+j)xk(x)}} & {e^{ - 2jxk(x)}} & {(1-j)e^{(1-j)xk(x)}} & { - 2jxk(x) - 3} \\ \end{matrix}\right] \end{aligned} Here, it is possible to check the determinant of ${\rm {\bf Q}}_{x_1 \to x_2 }$ given by $$\label{eq69} \left| {{\rm {\bf Q}}_{x_1 \to x_2 } } \right| = \exp \left( {\mathop{\rm tr}\left\{ {{\rm {\bf M}}_{x_1 \to x_2 } } \right\}} \right) = \exp \Big( { - 6\int_{x_1 }^{x_2 } {\frac{{k}'(x)}{k(x)}dx} } \Big) = \frac{k^6(x_1)}{k^6(x_2)},$$ in justification of (\ref{eq40}). To show the applicability of the approach, we take $a{ }_0(x) = - a^4x^{ - 4}$ for which (\ref{eq67}) becomes an Euler-Cauchy equation and therefore has an exact solution given by $$\label{eq70} f(x) = c_1 x^{m_1 } + c_2 x^{m_2 } + c_3 x^{m_3 } + c_4 x^{m_4 },$$ where $m_i = \textstyle{3 \over 2}\pm \textstyle{1 \over 2}\sqrt {5\pm 4\sqrt {1 + a^4} } ,i = 1,\ldots, 4$, and $c_i ,i = 1,\ldots, 4$ are arbitrary constants. Since $k(x) = ax^{ - 1}$, the kernel matrix (\ref{eq68}) simplifies into \label{eq71} \begin{aligned} &{\rm {\bf U}}(x) \\ &= \frac{ - 1}{2x}\left[ \begin{matrix} {2a - 3} & {(1+j)e^{(1-j)a}} & {e^{2a}} & {(1-j)e^{( {1 +j})a}} \\ {(1-j)e^{(j-1)a}} & {2ja - 3} & {(1+j)e^{(1+j)a}} & {e^{2ja}} \\ {e^{ - 2a}} & {(1-j)e^{ - (1+j)a}} & { - 2a - 3} & {(1+j)e^{( {j- 1})a}} \\ {(1+j)e^{ - (1+j)a}} & {e^{ -2ja}} & {(1-j)e^{(1-j)a}} &{ - 2aj - 3} \\ \end{matrix} \right] \\ &\equiv \frac{3}{2x}{\rm {\bf I}} + \frac{1}{x}{\rm {\bf N}} \\ \end{aligned}, where \textbf{N} is a traceless constant matrix with eigenvalues $\lambda _i = \pm \textstyle{1 \over 2}\sqrt {5\pm 4\sqrt {1 + a^4} }$, $i = 1,\ldots, 4$. For this case, the kernel matrix ${\bf U}(x)$ commutes with the transfer exponent ${\bf M}_{1\to x}$ so that (\ref{eq313}) becomes exact. By integration of (\ref{eq71}) we find $$\label{eq72} {\rm {\bf M}}_{1 \to x} = \frac{3\ln x}{2}\,{\rm {\bf I}} + \ln x\,{\rm {\bf N}}.$$ Obviously \textbf{I} and \textbf{N} commute, so that $$\label{eq73} {\rm {\bf Q}}_{1 \to x} = \exp \left( {{\rm {\bf M}}_{1 \to x} } \right) = x^{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-\nulldelimiterspace} 2}\exp \left( {\rm {\bf N}} \right) = x^{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-\nulldelimiterspace} 2}{\rm {\bf R}}\left[ {\exp \left( {\lambda _i \ln x} \right)\delta _{ij} } \right] {\rm {\bf R}}^{ - 1},$$ where the constant matrix \textbf{R} is the diagonalizer of \textbf{N}. Therefore, $\left| {{\rm {\bf Q}}_{1 \to x} } \right| = x^6$, as required by (\ref{eq69}). Furthermore, the elements of ${\rm {\bf Q}}_{1 \to x}$ involve linear combinations of $x^{r_i },i = 1,\ldots, 4$, for which $r_i = m_i$ hold. Finally from (\ref{eq6}), (\ref{eq7}), and (\ref{eq31}) we have $$\label{eq74} f(x) = \exp \left[ {{\rm {\bf \Phi }}(x)} \right]\,^t{\rm {\bf Q}}_{1 \to x} {\rm {\bf F}}\left( 1 \right),$$ where ${\rm {\bf \Phi }}(x)$ has become a constant vector, and ${\rm {\bf F}}\left( 1 \right)$ is a vector of arbitrary constants, to be determined by boundary conditions. However, comparing to (\ref{eq70}) shows that (\ref{eq74}) must be indeed the general solution of (\ref{eq67}). \subsection{Abel-Liouville-Ostogradski Formula} Found in 1827 by Abel for second-order differential equations and by Liouville and Ostogradsky in 1838 for the general case, the Wronskian should satisfy [26-28] $$\label{eq75} {W}'(x) + a_{n - 1} (x)W(x) = 0,$$ which implies that the Wronskian must be essentially a constant if $a_{n - 1} (x) \equiv 0$. Now we show that this equation can be readily reconstructed by DTMM. Following the discussions in sub-section 4.4 the Wronskian determinant takes the simple form $W = \left| {\rm {\bf D}} \right|\,\left| {\exp \left( {x{\rm {\bf K}}} \right)} \right|\,\,\left| {{\rm {\bf Q}}_{c \to x} } \right|\,\,\left| {\rm {\bf F}} \right|$, where $c \in \mathbb{S}$ is a constant and \textbf{F} is the matrix of independent vectors as defined in section 4.4. Therefore using (\ref{eq17}) and (\ref{eq39}) we find \label{eq78} \begin{aligned} W &= \prod_{i > j} {\left[ {k_i (x) - k_j (x)} \right]} \exp \left[ {x\mathop{\rm tr}\left\{ {\rm {\bf K}} \right\}} \right]\\ &\quad\times \exp \Big[ {ca_{n - 1} \left( c \right) - xa_{n - 1} (x) - \int_c^x {a_{n - 1} (x)dx} } \big]\prod_{i > j} {\frac{k_i \left( c \right) - k_j ( c)}{k_i (x) - k_j (x)}} \, \left| {\rm {\bf F}} \right| \\ &= \exp \left[ {ca_{n - 1} \left( c \right)} \right]\prod_{i > j} {\left[ {k_i \left( c \right) - k_j \left( c \right)} \right]} \,\left| {\rm {\bf F}} \right|\exp \Big[ { - \int_c^x {a_{n - 1} (x)dx} } \Big] \\ &\equiv A\exp \Big[ { - \int_c^x {a_{n -1}(x)dx} } \Big], \end{aligned} in which $A$ is a constant. Clearly the Wronskian given by (\ref{eq78}) satisfies the differential equation (\ref{eq75}). \section{Conclusions} We have presented a new analytical formalism for solution of linear homogeneous ordinary differential equations with variable coefficients. The formalism is based on the definition of jump transfer matrices and their extension into differential form. We presented the derivative lemma and the fundamental theorem of differential transfer matrix method, which support the exactness of the formalism. We also discussed the linear independency of solutions. The method is completely general, but fails in the presence of identical roots of characteristic equation. We have discussed a method to deal with corresponding singularities. The main advantage of the presented method is that it deals directly with the evolution of linear independent solutions, rather than their derivatives. The DTMM, as described in $\S 5$, when applied to wave propagation problems, leads to a novel approach for understanding the behavior of physical systems. \subsection*{Acknowledgement} One of the authors (Sina Khorasani) wishes to thank Prof. Beresford Parlett and Dr. Michael Parks at University of California at Berkeley for valuable discussions of this work. %\appendix \section{Appendix} %: Proof of (\ref{eq36})} We prove the validity of (\ref{eq36}), by showing first that both sides have the same derivative. From [41] the left-hand side of (\ref{eq36}) should obey the differential equation \label{eqA1} \begin{aligned} &\frac{d}{dx_2}\Big| \exp \Big[ {\int_{x_1}^{x_2} {{\rm {\bf H}}(t)^{ - 1}{\rm {\bf {H}'}}( t)dt} } \Big] \Big|\\ &= \Big| \exp \Big[{\int_{x_1}^{x_2} {{\rm {\bf H}}(t)^{ - 1}{\rm {\bf {H}'}}(t)dt} } \Big] \Big| \mathop{\rm tr}\left\{ {\rm {\bf H}}^{ - 1}\left( x_2\right){\rm {\bf H}}'\left( x_2 \right) \right\}. \end{aligned} Similarly, application of the derivative theorem for determinants [32, p.178] to the right-hand-side of (\ref{eq36}) results in $$\label{eqA2} \frac{d}{dx_2}\left| {\rm {\bf H}}^{ - 1}\left( x_1 \right){\rm {\bf H}}\left( x_2 \right) \right|=\frac{1}{\left|{\rm {\bf H}}\left( x_1 \right)\right|}\frac{d\left| {\rm {\bf H}}\left( x_2 \right) \right|}{dx_2}=\\ \frac{\left| {\rm {\bf H}}\left( x_2 \right) \right|}{\left|{\rm {\bf H}}\left(x_1 \right)\right|}\mathop{\rm tr}\left\{{\bf H}^{-1}(x_2){\bf H}'(x_2)\right\}.$$ Therefore, both sides satisfy the same differential equation, and thus $$\label{eqA3} \Big| \exp \Big[ {\int_{x_1}^{x_2} {{\rm {\bf H}}(x)^{ - 1}{\rm {\bf {H}'}}(x)dx} } \Big] \Big| = w(x_1)\left| {\rm {\bf H}}^{ - 1}\left( x_1 \right){\rm {\bf H}}\left( x_2 \right) \right|,$$ where $w(\cdot)$ is a function of only $x_1$. But upon setting $x_1=x_2$ in (\ref{eqA3}), one gets $w(\cdot)\equiv1$ and hence the claim. \subsection*{Matrix Exponential of 2$\times$2 Matrices} The difficulties in evaluation of matrix exponential have been well known for a long time [54,55]. While computation of matrix exponential is possible by general [56-58] and approximate [59] methods, simple analytical expression exists only for few cases including the Euler-Rodrigues formula [59] for $3\times 3$ skew-symmetric matrices. Here, we report an exact expansion of matrix exponential for arbitrary $2\times 2$ square matrices. Assume that \textbf{M} is a square matrix with arbitrary complex elements. We define the traceless matrix $\textbf{A}=\textbf{M}- \textstyle{1 \over 2}\mathop{\rm tr}\left\{ {\rm {\bf M}} \right\}{\rm {\bf I}}$. Since \textbf{A} and $\textstyle{1 \over 2} \mathop{\rm tr}\left\{ {\rm {\bf M}} \right\}{\rm {\bf I}}$ obviously commute, one can deduce that $$\label{eq76} \exp \left( {\rm {\bf M}} \right) = \exp \left( {\textstyle{1 \over 2}\mathop{\rm tr}\left\{ {\rm {\bf M}} \right\}} \right)\exp \left( {\rm {\bf A}} \right).$$ Now, \textbf{A} satisfies the property ${\rm {\bf A}}^2 = \left| {\rm {\bf A}} \right|\,{\rm {\bf I}}$, so that [60] $$\label{eq77} \exp \left( {\rm {\bf A}} \right) = \cos \left( \delta \right)\;{\rm {\bf I}} + {\rm sinc}\left( \delta \right)\;{\rm {\bf A}},$$ in which $\delta = \sqrt { - \left| {\rm {\bf A}} \right|}$ and $\mathop{\rm sinc}(x)=\sin(x)/x$. Finally, $$\label{eq79} \exp \left( {\rm {\bf M}} \right) = \exp \left( \textstyle{1 \over 2}\mathop{\rm tr}\left\{ {\rm {\bf M}} \right\} \right)\,\left[ {\cos (\delta) \;{\rm {\bf I}} + \mathop{\rm sinc}(\delta) \;{\rm {\bf A}}} \right].$$ \begin{thebibliography}{00} \bibitem{1} M. Braun, \textit{Differential Equations and Their Applications}, 4th ed., Springer-Verlag, New York, 1993. \bibitem{2} N. H. Ibragimov, \textit{Elementary Lie Group Analysis and Ordinary Differential Equations}, John Wiley {\&} Sons, Chichester, 1999. \bibitem{3} O. Gat, \textit{Symmetry algebras of third-order ordinary differential equations}, J. Math. Phys., 33 (1992), pp. 2966-2971. \bibitem{4} L. M. Berkovich and S. Y. Popov, \textit{Group analysis of ordinary differential equations of the order of n}$>$\textit{2}, Sym. Nonlin. Math. Phys., 1 (1997), pp. 164-171. \bibitem{5} W. R. Oudshoorn and M. Van Der Put, \textit{Lie symmetries and differential Galois groups of linear equations}, Math. Comp., 71 (2001), pp. 349-361. \bibitem{6} L. M. Berkovich, \textit{The method of an exact linearization of n-order ordinary differential equations}, J. Nonlin. Math. Phys., 3 (1996), pp. 341-350. \bibitem{7} R. M. Edelstein, K. S. Govinder and F. M. Mahomed, \textit{Solutions of ordinary differential equations via nonlocal transformations}, J. Phys. A: Math. Gen., 34 (2001), pp. 1141-1152. \bibitem{8} L. M. Berkovich, \textit{Transformations of ordinary differential equations: local and nonlocal symmetries}, Proc. Inst. Math. NAS Ukraine, 30 (2000), pp. 25-34. \bibitem{9} M. Bronstein and S Lafaille, \textit{Solutions of linear ordinary differential equations in terms of special functions}, in Proceedings of ISSAC'2002, ACM Press, Lille, 2002, pp. 23-28. \bibitem{10} A. Holubec and A. D. Stauffer, \textit{Efficient solution of differential equations by analytic continuation}, J. Phys. A: Math. Gen., 18 (1985), pp. 2141-2149. \bibitem{11} E. A. Coutsias, T. Hagstorm and D. Torres, \textit{An efficient spectral method for ordinary differential equations with rational function coefficients}, Math. Comp., 65 (1996), pp. 611-635. \bibitem{12} E. Babolian and M. M. Hosseini, \textit{A modified spectral method for numerical solution of ordinary differential equations with non-analytic solution}, Appl. Math. Comp., 132 (2002), pp. 341-351. \bibitem{13} W. N. Everitt, K. H. Kwon, L. L. Littlejohn and R. Wellman, \textit{Orthogonal polynomial solutions of linear ordinary differential equations}, J. Comp. Appl. Math., 133 (2001), pp. 85-109. \bibitem{14} V. Kozlov and V. Maz'ya, \textit{Differential Equations with Operator Coefficients}, Springer, Berlin, 1999. \bibitem{15} V. K. Dzyadyk, \textit{Approximation Methods for Solutions of Differential and Integral Equations}, VSP, Utrecht, 1995. \bibitem{16} B. Haile, \textit{Analytic solution of n-th order differential equations at a singular point}, Elec. J. Diff. Eq., 2002 (2002), pp. 1-14. \bibitem{17} W. T. Reid, \textit{Ordinary Differential Equations}, John Wiley {\&} Sons, New York, 1971, p. 58. \bibitem{18} S. Khorasani and K. Mehrany, \textit{Analytical solution of wave equation for arbitrary non-homogeneous media}, Proc. SPIE, 4772 (2002), pp. 25-36. \bibitem{19} S. Khorasani and K. Mehrany, \textit{Differential transfer matrix method for solution of one-dimensional linear non-homogeneous optical structures}, J. Opt. Soc. Am. B, 20 (2003), pp. 91-96. \bibitem{20} P. Yeh, \textit{Optical Waves in Layered Media}, Wiley, New York, 1988. \bibitem{21} A. Fassler and E. Stiefel, \textit{Group Theoretical Methods and Their Applications}, Birkhauser, Boston, 1992. \bibitem{22} K. Mehrany and S. Khorasani, \textit{Analytical solution of non-homogeneous anisotropic wave equation based on differential transfer matrices}, J. Opt. A: Pure Appl. Opt., 4 (2002), pp. 524-635. \bibitem{23} D. Sarafyan, \textit{Approximate solution of ordinary differential equations and their systems through discrete and continuous embedded Runge-Kutta formulae and upgrading of their order}, Comp. Math. Appl., 28 (1992), pp. 353-384. \bibitem{24} J. C. Butcher, \textit{Numerical methods for ordinary differential equations in the 20th century}, J. Comp. Appl. Math., 125 (2000), pp. 1-29. \bibitem{25} W. Kahan and R. Li, \textit{Composition constants for rising the orders of unconventional schemes for ordinary differential equations}, Math. Comp., 66 (1997), pp. 1089-1099. \bibitem{26} S. K. Godunov, \textit{Ordinary Differential Equations with Constant Coefficients}, American Mathematical Society, Providence, 1997. \bibitem{27} I. G. Petrovski, \textit{Ordinary Differential Equations}, Prentice-Hall, Englewood Cliffs, 1966, p. 109. \bibitem{28} T. Apostol, \textit{Calculus}, 2nd ed., Blaisdell, Waltham, 1969, vol. 2, p. 161. \bibitem{29} E. G. C. Poole, \textit{Introduction to the Theory of Linear Differential Equations}, Clarendon Press, Oxford, 1936, p. 12. \bibitem{30} P. Lancaster, \textit{Theory of Matrices}, Academic Press, New York, 1969, p. 36. \bibitem{31} S. Khorasani and A. Adibi, \textit{New analytical approach for computation of band structure in one-dimensional periodic media}, Opt. Commun., 216 (2003), pp. 439-451. \bibitem{32} J. R. Mangus, and H. Neudecker, \textit{Matrix Differential Calculus with Applications in Statistics and Econometrics}, John Wiley {\&} Sons, Chichester, 1988, p. 151. \bibitem{33} F. J. Dyson, \textit{The S matrix in quantum electrodynamics}, Phys. Rev., 75 (1949) pp. 1736-1755. \bibitem{34} P. Roman, \textit{Advanced Quantum Theory}, Addison-Wesley, Palo Alto, 1965, p. 310. \bibitem{35} J. D. Bjorken, and S. D. Drell, \textit{Relativistic Quantum Mechanics}, McGraw-Hill, 1965, p. 177. \bibitem{36} S. Khorasani, \textit{Reply to Comment on `Analytical solution of non-homogeneous anisotropic wave equation based on differential transfer matrices"'}, J. Opt. A: Pure Appl. Opt., 5 (2003), pp. 434-435. \bibitem{37} S. Khorasani and A. Adibi, \textit{New Matrix Method for Analytical Solution of Linear Ordinary Differential Equations}, to be published in Proc. 4th Dynamic Systems and Applications, Dynamic Publishers, Atlanta, 2003. \bibitem{38} J.-C. Evard, \textit{On matrix functions which commute with their derivative}, Linear Algebra Appl., 68 (1985), pp. 145-178. \bibitem{39} J.-C. Evard, \textit{Invariance and commutativity properties of some classes of solutions of the matrix differential equation $X(t)X'(t)=X'(t)X(t)$}, Linear Algebra Appl., 218 (1995), pp. 89-102. \bibitem{40} N. P. Erugin, \textit{Linear Systems of Ordinary Differential Equations}, Academic Press, New York, 1966, p. 38. \bibitem{41} J. Cornin, \textit{Differential Equations}, 2nd ed., Marcel Dekker, New York, 1994, p. 73. \bibitem{42} O. Cormier, M. F. Singer, B. M. Trager, and F. Ulmer, \textit{Linear differential operators for polynomial equations}, J. Sym. Comp., 34 (2002) pp. 355-398. \bibitem{43} V. Neagoe, \textit{Inversion of the Van der Monde matrix}, IEEE Sig. Proc. Lett., 3 (1996), pp. 119-120. \bibitem{44} J. C. F. Sturm, \textit{M$\acute e$moire sur les $\acute e$quations differentieles lin$\acute e$aire de second ordre}, J. Math. Pures Appl., 1 (1836), pp. 106-186. \bibitem{45} C. A. Swanson, \textit{Comparison and Oscillation Theory of Linear Differential Equations}, Academic Press, New York, 1968. \bibitem{46} W. Coppel, \textit{Disconjugacy}, Springer-Verlag, New York, 1971. \bibitem{47} P. Hartman, \textit{Ordinary Differential Equations}, 2nd ed., Birkhouse, Boston, 1982. \bibitem{48} Z. Deng, and S. Ruan, \textit{On second order linear differential systems}, in \textit{Differential Equations and Control Theory}, edited by Z. Deng, Z. Liang, G. Lu and S. Ruan, Marcel Dekker, New York, 1996, p. 21. \bibitem{49} P. W. Eloe and J. Henderson, \textit{Positive solutions for higher order ordinary differential equations}, Elec. J. Diff. Eq., 1995 (1995), pp. 1-8. \bibitem{50} M. Naito, \textit{On the number of zeros of nonoscillatory solutions to higher-order linear ordinary differential equations}, Monatsh. Math., 136 (2002), pp. 237-242. \bibitem{51} T. Regge and J. A. Wheeler, \textit{Stability of a Schwarzschild singularity}, Phys. Rev., 108 (1957) pp. 1063-1069. \bibitem{52} S. Mano, H. Suzuki, and E. Takasugi, \textit{Analytic solutions of the Regge-Wheeler equation and the post-Minkowskian expansion}, Prog. Theor. Phys., 96 (1996) pp. 549-565. \bibitem{53} S. Khorasani, and B. Rashidian, \textit{Modified transfer matrix method for conducting interfaces}, J. Opt. A: Pure Appl. Opt., 4 (2002) pp. 251-256. \bibitem{54} C. Moller and C. F. Van Loan, \textit{Nineteen dubious ways to compute the exponential of a matrix}, SIAM Rev., 20 (1978), pp. 801-836. \bibitem{55} M. Parks, \textit{A Study of Algorithms to Compute the Matrix Exponential}, PhD Dissertation, Department of Applied Mathematics, University of California at Berkeley, 1994. \bibitem{56} I. E. Leonard, \textit{The matrix exponential}, SIAM Rev., 38 (1996), pp. 507-512. \bibitem{57} E. Liz, \textit{A note on the matrix exponential}, SIAM Rev., 40 (1998), pp. 700-702. \bibitem{58} W. A. Harris, J. P. Fillmore, and D. R. Smith, \textit{Matrix exponential-another approach}, SIAM Rev., 43 (2001), pp. 694-706. \bibitem{59} A. Zanna and H. Z. Munthe-Kaas, \textit{Generalized polar decompositions for the approximation of the matrix exponential}, SIAM J. Matrix Anal. Appl., 23 (2002), pp. 840-862. \bibitem{60} I. M. Mladenov, \textit{Saxon-Hunter theorem via matrix exponential}, Europhys. Lett., 37 (1997), pp. 159-164. \end{thebibliography} \end{document}