Initial conditions for second order ODE with complex stiffness

Question

I tried this on Math Stack Exchange. I'm trying to find initial conditions to ensure systems of the form stay bounded

$$\ddot{x}_i+\sum_{j=1}^N k_{ij} x_j = 0, \quad k_{ij} \in \mathbb{C}.$$

For simplicity let's say the $k_{ij}$ lie in the first quadrant. I've managed to figure it out for 1D

$$\ddot{x}+kx=0 \quad k \in \mathbb{C}.$$

This equation has a solution of the form

$$x(t)= a e^{\sqrt{-k}t} + b e^{-\sqrt{-k}t}.$$

Since $k$ is complex the first term diverges exponentially. The initial condition

$$x'(0)+\sqrt{-k} x(0)=0$$

kills off the first term, making the solution finite for all $t>0$.

I'm trying to generalize this idea to a set of coupled oscillators. If we assume a solution of the form

$$x_i(t)=r_{ik} e^{\lambda_k t},$$

and insert it into the top equation, we find the $-\lambda_k^2$ are the eigenvalues of the matrix $\textbf{k}$, which has elements $k_{ij}$, and they have associated right eigenvectors with elements $r_{ik}$. The $\lambda_k$ are given by solutions to the characteristic equation

$$|\textbf{k} + \lambda_k^2 \textbf{I}|=0.$$

Thus, the most general solution to the top equation can be written as

$$x_i(t)=\sum_{k=1}^N r_{ik} (a_k e^{\lambda_k t} + b_k e^{-\lambda_k t} ) $$

So we have to find a set of initial constraints that ensures the $a_k=0$. Taking the derivative, and setting t=0, we find the two equations

$$x_i(0)=\sum_{k=1}^N r_{ik} (a_k + b_k )$$

$$x_i'(0)=\sum_{k=1}^N r_{ik} \lambda_k (a_k - b_k )$$

We can solve for the quantities $a_k + b_k$ and $a_k - b_k$ by finding the inverse of the eigenvector matrix $\textbf{r}$, which has elements $r_{ik}$. Multiplying the two initial condition equations by $\textbf{r}^{-1}$ we get

$$\sum_{i=1}^N r_{ji}^{-1} x_i(0)=a_j+b_j$$

$$\frac{1}{\lambda_j}\sum_{i=1}^N r_{ji}^{-1} x_i'(0)=a_j-b_j$$

Adding these two equations together cancels out the $b_j$. Setting the $a_j=0$ we find our set of initial conditions

$$\sum_{i=1}^N r_{ji}^{-1} x_i'(0)+\lambda_j \sum_{i=1}^N r_{ji}^{-1} x_i(0) = 0 $$

Now all this seems alright to me, but maybe I made a mistake somewhere. I'm just a simple physicist, and sometimes I get confused when summing over multiple indices. I wrote some Mathematica code to try to implement all this, and it does not produce the desired result. If you run it you will see it does not kill off the diverging terms. It does however work quite nicely if you set N=1 (I called it nm in the program). Is my derivation wrong, or is my code wrong, or both?

nm = 2;
x = ToExpression["x" <> ToString[#]][t] & /@ Range[nm];
k = RandomComplex[2 + 2*I, {nm, nm}];
eqns = D[x[[#]], t, t] + Sum[k[[#, i]]*x[[i]], {i, 1, nm}] == 0 & /@ 
   Range[nm];
eigVals = Eigenvalues[k];
\[Lambda] = Sqrt[-eigVals];
eigVecs = Eigenvectors[k];
eigVecsInverse = Inverse[eigVecs];
initsa = 0 == (Sum[
        eigVecsInverse[[#, i]]*D[x[[i]], t] /. t -> 0, {i, 1, 
         nm}] + \[Lambda][[#]]*
        Sum[eigVecsInverse[[#, i]]*x[[i]] /. t -> 0, {i, 1, nm}]) & /@
    Range[nm];
initsb = Thread[(x /. t -> 0) == 1 + I];
sys = Join[eqns, initsa, initsb];
aSol = DSolve[sys, x, t];
TrigToExp[aSol]