We extend the Bloch-decomposition based time-splitting spectral method introduced in an earlier paper [Z. Huang, S. Jin, P. Markowich, C. Sparber, A Bloch decomposition based split-step pseudo spectral method for quantum dynamics with periodic potentials, SIAM J. Sci. Comput. 29 (2007) 515–538] to the case of (non-)linear Klein–Gordon equations. This provides us with an unconditionally stable numerical method which achieves spectral convergence in space, even in the case where the periodic coefficients are highly oscillatory and/or discontinuous. A comparison to a traditional pseudo-spectral method and to a finite difference/volume scheme shows the superiority of our method. We further estimate the stability of our scheme in the presence of random perturbations and give numerical evidence for the well-known phenomenon of Anderson’s localization.