A numerical algorithm is presented that simulates large deformations of heterogeneous, viscoelastic materials in two dimensions. The algorithm is based on a spectral/finite difference method and uses the Eulerian formulation including objective derivatives ofthe stress tensor in the rheological equations. The viscoelastic rheology is described bythe linear Maxwell model, which consists of an elastic and viscous element connected inseries. The algorithm is especially suitable to simulate periodic instabilities. The derivatives in the direction of periodicity are approximated by spectral expansions, whereas the derivatives in the direction orthogonal to the periodicity are approximated by finite differences. The 1‐D Eulerian finite difference grid consists of centre and nodal points and has variable grid spacing. Time derivatives are approximated with finite differences using an implicit strategy with a variable time step. The performance of the numerical code is demonstrated by calculation, for the first time, of the pressure field evolution during folding of viscoelastic multilayers. The algorithm is stable for viscosity contrasts up to 5 × 105, which demonstrates that spectral methods can be used to simulate dynamical systems involving large material heterogeneities. The successful simulations show that combined spectral/finite difference methods using the Eulerian formulation are a promising tool to simulate mechanical processes that involve large deformations, viscoelastic rheologies and strong material heterogeneitie
