Given a discrete-time linear time-invariant (LTI) system like
$$ x(k+1) = A x(k) $$
the solution at a given instant $k$ can be computed using the state transition matrix $A^k$ by
$$ x(k) = A^k x(0) $$
For continuous-time LTI systems, the state transition matrix can be computed using the matrix exponential, which gives the continuous state transition matrix with high accuracy.
However, for discrete-time systems, the state transition matrix is not computed with the matrix exponential but with $A^k$. In Matlab, for example, I could compute $A^k$ directly using a syntax like
A^k, which I guess just performs $k$ subsequent matrix multiplications. However, I am not sure if this is the best approach if $k$ is large.
Question: What is the most efficient way (in terms of accuracy) to compute $A^k$ for large $k$ numerically (preferable in Matlab)?