We appreciate S. M. Shahruz for his interest in our paper and for the useful comments on the numerical method for solving the transformed system
$Y¨t+eCtΛ2−C2e−CtYt=0,$

$Y0=X0,Y˙0=CX0+X˙0$
(1)
As pointed out by Shahruz, the numerical method proposed in our paper is an approximate approach. The approximation comes from the assumption that the time-varying matrix $eCtΛ2−C2eCt$ is time-invariant within short time duration $jh$j=0,1,2,…,N$ where $h$ is the time step, i.e.,
$eCt[Λ2−C2]e−Ct≈ejhC[Λ...$