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)?
Some general comments:
The matrix $A^k$ is not what most people refer to as the state transition matrix. Instead this is the matrix $A$, since it transfers the state at time $k$ to the state at time $k+1$. In general, for the discrete linear system $x_{k+1} = A_kx_{k}$ you can write the matrix $$\Phi(k,k') = \prod_{j=k'}^{k-1}A_j$$ so that $$x_k = \Phi(k,k')x_{k'}.$$ Note that in the time-invariant case we have $\Phi(k,k') = A^{k-k'}$ and $\Phi(k,0) = A^k$, which is what you have. Sometimes $\Phi(k,k')$ is referred to as the state transition matrix and $A$ is referred to as its generator.
Programming:
The best way to have matlab simulate the solution to $x_{k+1} = A_kx_k$, is to use the equation itself rather than trying to calculate any transition matrix (other than for a single step). You can do this by storing the values of the state vector and computing the next value using the update equation. If you really need an analytic formula for $\Phi(k,k')$ or $A^k$, you can find the Jordan form and take powers of the Jordan blocks, which have well-known formulas that depend on the block type.