We derive a numerical scheme to compute invariant manifolds for timevariant discrete dynamical systems, i.e., nonautonomous difference equations. Our universally applicable method is based on a truncated Lyapunov-Perron operator and computes invariant manifolds using a system of nonlinear algebraic equations which can be solved both locally using (nonsmooth) inexact Newton, and globally using continuation algorithms. Compared to other algorithms, our approach is quite flexible, since it captures time-dependent, nonsmooth, noninvertible or implicit equations and enables us to tackle the full hierarchy of strongly stable, stable and center-stable manifolds, as well as their unstable counterparts.Our results are illustrated using a test example and are applied to a population dynamical model and the Hénon map. Finally, we discuss a linearly implicit Euler-Bubnov-Galerkin discretization of a reaction diffusion equation in order to approximate its inertial manifold.