The wrinkling phenomenon of a compressed elastic thin film bonded to a viscous layer is studied. Linear stability analysis (LSA) shows that, for materials with cubic crystalline symmetry, the sign of the degree of elastic anisotropy ζ = C 12 +2C 44 C 11 −1 plays an important role in the anisotropy of the buckling instability of the thin film system. More precisely, the growth rate of the fastest growing wave number, taking as a function of directions, reaches a peak in the < 100 > directions for ζ > 0, and in the < 110 > directions for ζ < 0. A highly efficient semi-implicit spectral method is established. The numerical experiments of long time wrinkling evolution processes of a 1+2 dimensional system verified the LSA results, successfully simulated anisotropic wrinkling pattern formation and coarsening, produced a power law scaling and reproduced certain featured phenomena observed in physical experiments.