Variational integrators provide a systematic way to derive geometric numerical methods for Lagrangian dynamical systems, which preserve a discrete multisymplectic form as well as momenta associated to symmetries of the Lagrangian via Noether's theorem. An inevitable prerequisite for the derivation of variational integrators is the existence of a variational formulation for the considered dynamical system. Even though this is the case for a large class of systems, there are many interesting examples which do not belong to this class, e.g., equations of advection-diffusion type like they are often found in fluid dynamics or plasma physics. We propose the application of the variational integrator method to so called formal Lagrangians, which allow us to embed any dynamical system into a Lagrangian system by doubling the number of variables. Thereby we are able to derive variational integrators for arbitrary systems, extending the applicability of the method significantly. A discrete version of the Noether theorem for formal Lagrangians yields the discrete momenta preserved by the resulting numerical schemes. The theory is applied to dynamical systems from fluid dynamics and plasma physics like the vorticity equation, the Vlasov-Poisson system and magnetohydrodynamics, including numerical examples. Recent attempts of applying this method also to noncanonical Hamiltonian ODEs will be sketched.