Let's find directly the asymptotic result of Prop.1.
Diagonalize the matrix $A_0$, with $D$ the diagonal matrix having diagonal elements equal to $(\lambda_i)$ for $i=1,...,d$ and $U$ the matrix of eigenvectors.
$$A_0 = U^{-1}DU$$
The system of equation is then equivalent to
$$\begin{align}
&\Longleftrightarrow \frac{d }{d t}\mathbf{y} = \left(U^{-1}DU + A_1 t^{-1} + A_2 t^{-2} + \dots A_n t^{-n} \right) \mathbf{y}\\
&\Longleftrightarrow \frac{d }{d t}(U\mathbf{y}) = \left(DU + UA_1 t^{-1} + UA_2 t^{-2} + \dots UA_n t^{-n} \right) \mathbf{y} \\
&\Longleftrightarrow \frac{d }{d t}(U\mathbf{y}) = \left(D + UA_1U^{-1} t^{-1} + UA_2U^{-1} t^{-2} + \dots UA_nU^{-1} t^{-n} \right) (U\mathbf{y}) \tag{1}\\
\end{align}$$
Let's denote $\mathbf{z} = U\mathbf{y}$, then
$$(1) \Longleftrightarrow \frac{d }{d t}\mathbf{z} = \left(D + UA_1U^{-1} t^{-1} + UA_2U^{-1} t^{-2} + \dots UA_nU^{-1} t^{-n} \right) \mathbf{z} \tag{2}$$
The $i^{\text{th}}$ equation of $(2)$, for $i=1,...,d$ has the following form
$$\frac{d z_i }{d t}(t) = \lambda_i z_i(t)+ t^{-1}\sum_{j=1}^d (\alpha_{1j}z_{j}(t)) + \color{red}{t^{-2}}\sum_{j=1}^d (\alpha_{2j}z_{j}(t)) + \dots + \color{red}{t^{-n}}\sum_{j=1}^d (\alpha_{nj}z_{i}(t)) \tag{3} $$
For $t$ large ($t\to +\infty$), the terms in red color become small. By the same argument as for the unidimensional case, we can ignore these terms, and then, for $i=1,...,d$, we have
$$(3) \implies \frac{d z_i }{d t}(t) \sim \lambda_i z_i(t) + \frac{1}{t}\sum_{j=1}^d (\alpha_{1j}z_{j}(t)) \qquad \forall i=1,...,d \tag{4}$$
Let's solve the system of equations $(4)$, we have
$$\begin{align}
(4) &\Longleftrightarrow \left(e^{-\lambda_i t} z_i(t) \right)' = \frac{1}{t}\sum_{j=1}^d \left(\alpha_{ij} e^{-\lambda_i t}z_{j}(t) \right) \qquad \forall i=1,...,d
\end{align}$$
For the sake of simplicity, let's denote $p_i(t) = e^{-\lambda_i t} z_i(t)$, and $ \mathbf{p} =\left(p_1(t),...,p_d(t) \right)$. The matrix form of the relation between $p_i(t)$ and $z_i(t)$ is
$$\mathbf{p} = \Gamma^{-1}\mathbf{z}$$
with $\Gamma$ - a diagonal matrix having elements $e^{\lambda_i t}$ for $i=1,...,d$
And the matrix form of the new system of equation is
$$\frac{d}{dt} \mathbf{p} = \frac{1}{t} UA_1U^{-1} \mathbf{p} \tag{5}$$
Diagonalize the matrix $A_1$ with $H$ the diagonal matrix having diagonal elements $\beta_1,...,\beta_d$ and $V$ the matrix of eigenvectors.
$$A_1 = V^{-1}HV$$
then
$$\begin{align}
(5) &\Longleftrightarrow \frac{d}{dt}\mathbf{p} = \frac{1}{t} UV^{-1}HVU^{-1} \mathbf{p} \\
&\Longleftrightarrow \frac{d}{dt}(VU^{-1}\mathbf{p} )= \frac{1}{t} H (VU^{-1} \mathbf{p}) \tag{6}\\
\end{align}$$
$(6)$ can be solved easily. Indeed, each equation of $(6)$ has the following form
$$\frac{d}{dt}q_i(t) = \frac{\beta_i}{t}q_i(t) \Longleftrightarrow q_i(t) = c_it^{\beta_i} \qquad \forall i=1,...,d $$
Then
$$\begin{align}
(6)&\Longleftrightarrow VU^{-1}\mathbf{p} = \begin{pmatrix}
c_1t^{\beta_1} \\
c_2t^{\beta_2} \\
\cdot \cdot \cdot\\
c_dt^{\beta_d} \\
\end{pmatrix} \\
&\Longleftrightarrow
VU^{-1}\Gamma^{-1}\mathbf{z} = \begin{pmatrix}
c_1t^{\beta_1} \\
c_2t^{\beta_2} \\
\cdot \cdot \cdot\\
c_dt^{\beta_d} \\
\end{pmatrix} \\
&\Longleftrightarrow
\Gamma^{-1}VU^{-1}\mathbf{z} = \begin{pmatrix}
c_1t^{\beta_1} \\
c_2t^{\beta_2} \\
\cdot \cdot \cdot\\
c_dt^{\beta_d} \\
\end{pmatrix} \qquad \text{as } \Gamma \text{ is diagonal matrix}\\
&\Longleftrightarrow
VU^{-1}U\mathbf{y} = \Gamma\begin{pmatrix}
c_1t^{\beta_1} \\
c_2t^{\beta_2} \\
\cdot \cdot \cdot\\
c_dt^{\beta_d} \\
\end{pmatrix} \\
&\Longleftrightarrow
\color{red}{\mathbf{y} = V^{-1} \begin{pmatrix}
c_1t^{\beta_1} e^{\lambda_1 t}\\
c_2t^{\beta_2}e^{\lambda_2 t} \\
\cdot \cdot \cdot\\
c_dt^{\beta_d}e^{\lambda_d t} \\
\end{pmatrix}} \tag{7} \\
\end{align}$$
with $V$ the matrix of eigenvectors defined by $A_1 = VHV^{-1}$.
$(7)$ is the more accurate solution than the one in your Prop.1. If we seek only the most major term for each $y_i(t)$, then we can do as follows
$$\color{red}{y_i(t) = \sum_{j=1}^d v_{ij}c_jt^{\beta_j}e^{\lambda_j t}} \sim v_{i\color{blue}{g}(i)}c_{\color{blue}{g}(i)} t^{\beta_{\color{blue}{g}(i)}}e^{\lambda{_\color{blue}{g}(i)} t} \qquad \forall i=1,...,d $$
with $\color{blue}{g}(i) =\underset{j=1,...,d}{\operatorname{argmax}} \{\lambda_j | v_{ij}c_j \ne 0\}$.
Remark: if all the coefficients $v_{ij}, c_j$ for $1\le i,j \le d$ are not $0$, then for all $i=1,...,d$, we have $\color{blue}{g}(i) = \color{blue}{g} =\underset{j=1,...,d}{\operatorname{argmax}} \{\lambda_j \}$ be the index $j^{th}$ such that $\lambda_j=\color{purple}{\lambda}$ (the P.E eigenvalue), and so we have
$$y_i(t) \sim v_{i\color{blue}{g}}c_\color{blue}{g} t^{\beta_\color{blue}{g}}e^{\color{purple}{\lambda} t} \qquad \forall i=1,...,d $$