  We use a modification of the parameterization method to study invariant
manifolds for difference equations. We establish existence, regularity, smooth
dependence on parameters and study several singular limits, even if the
difference equations do not define a dynamical system. This method also leads
to efficient algorithms that we present with their implementations. The
manifolds we consider include not only the classical strong stable and unstable
manifolds but also manifolds associated to non-resonant spaces.
  When the difference equations are the Euler-Lagrange equations of a discrete
variational we present sharper results. Note that, if the Legendre condition
fails, the Euler-Lagrange equations can not be treated as a dynamical system.
If the Legendre condition becomes singular, the dynamical system may be
singular while the difference equation remains regular.
  We present numerical applications to several examples in the physics
literature: the Frenkel-Kontorova model with long-range interactions and the
Heisenberg model of spin chains with a perturbation. We also present extensions
to finite differentiable difference equations.
