We present a lattice formulation of an interaction phi/Lambda F (F) over tilde between an axion and some U(1) gauge sector with the following properties: it reproduces the continuum theory up to O(dx(mu)(2)) corrections, it preserves exact gauge invariance and shift symmetry on the lattice, and it is suitable for self-consistent expansion of the Universe. The lattice equations of motion can no longer be solved by explicit methods, but we propose an implicit method to overcome this difficulty, which preserves the relevant system constraints down to arbitrary (tunable) precision. As a first application we study, in a comoving grid in (3 +1) dimensions, the last efolds of axion-inflation with quadratic potential and the preheating stage following afterwards. We fully account for the inhomogeneity and non-linearity of the system, including the gauge field contribution to the expansion rate of the Universe and its backreaction into the axion dynamics. We characterize in detail, as a function of the coupling, the energy transfer from the axion to the gauge field. Two coupling regimes are identified, sub- and super-critical, depending on whether the final energy fraction stored in the gauge field is below or above similar to 50% of the total energy. The Universe is very efficiently reheated for super-critical couplings, rapidly entering in a radiation dominated stage. Our results on preheating confirm previously published results.