Consider the linear system
$$\frac{du}{dt}=Ju, \qquad u(0) = u_0$$
where $J \in \mathbb R^{2 \times 2}$ has a pair of complex conjugate eigenvalues $a \pm i b$ and initial condition $u_0 \in \mathbb R^2$. I want to show that the solution obtained via the matrix exponential is real.
Suppose $J$ is diagonalizable, i.e., $J = T \Lambda T^{-1},$ then we have
$$u(t)=e^{Jt}u_0 = T e^{\Lambda t} T^{-1}u_0 = T \begin{bmatrix} e^{(a+bi)t} & 0 \\ 0 & e^{(a-bi)t} \end{bmatrix} T^{-1}u_0 = e^{at}T \begin{bmatrix} \cos(bt)+i\sin(bt) & 0 \\ 0 & \cos(bt)-i\sin(bt) \end{bmatrix} T^{-1}u_0 = e^{at}(\cos(bt)TIT^{-1}u_0 + \sin(bt)iT \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix} T^{-1}u_0) = e^{at} (\cos(bt)u_0 + \sin(bt) iT \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix}T^{-1}u_0)$$
$u_0, \cos(bt), \sin(bt), e^{at}$ are all reals, therefore, we want to show that $iT \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix}T^{-1}$ is real. But I'm not quite sure how to do that...
I don't know how to get that $u(t)$ is real following your procedure. Nonetheless, if $J \in \mathbb R^{2\times 2}$ (otherwise, no chance that you will have $u(t) \in \mathbb R^2$), you can write the series expansion of $\exp(Jt)$ and obtain $$ u(t) = e^{Jt}u_0 = \left(\sum_{k=0}^\infty \frac{J^kt^k}{k!}\right) u_0. $$ The series is guaranteed to converge in the space of matrices with real coefficients (see, e.g., here) for all $t > 0$. Therefore, $u(t)$ is real.
Addendum: you should have a look here.