The law of the iterated logarithm (LIL) for the time-homogeneous Markov process with a unique invariant measure characterizes the almost sure maximum possible fluctuation of time averages around the ergodic limit. Whether a numerical approximation can preserve this asymptotic pathwise behavior remains an open problem. In this work, we give a positive answer to this question and establish the LIL for the numerical approximation of such a process under verifiable assumptions. The Markov process is discretized by a decreasing time-step strategy, which yields the non-homogeneous numerical approximation but facilitates a martingale-based analysis. The key ingredient in proving the LIL for such numerical approximation lies in extracting a quasi-uniform time-grid subsequence from the original non-uniform time grids and establishing the LIL for a predominant martingale along it, while the remainder terms converge to zero. Finally, we illustrate that our results can be flexibly applied to numerical approximations of a broad class of stochastic systems, including SODEs and SPDEs.