We present LaCoDe (Lagrangian Compressible Deformation), a MATLAB solver for the Stokes equations for compressible non-Newtonian visco-elastic flow in two dimensions, based on a Lagrangian formulation of the Finite Element Method. The incompressible Boussinesq approximation is a widespread assumption in numerical models of lithospheric deformation, thus potentially masking a significant contribution of mechanisms linked to volumetric changes that occur in the asthenospheric mantle and the lithosphere. LaCoDe employs a compressible formulation of the Stokes equations designed to address such volume-changing processes. First, we provide a description of the equations governing the deformation of Earth rocks and detailed overview of the algorithm, its numerical implementation, treatment of the non-linearities rising from the compressible formulation and ineleastic deformation, and the remeshing algorithm that tracks and transfers the physical fields from a highly-distorted mesh to a high-quality one. LaCoDe is then benchmarked by comparing numerical results to analytical solutions for the bending of a thin elastic beam under a constant uniform load, flow around a rigid inclusion, Rayleigh-Taylor instability, stress build-up in a visco-elastic Maxwell body, and Couette flow with viscous heating. The Rayleigh-Taylor instability test is further used to demonstrate the accuracy of the remeshing algorithm. The importance of including volumetric strain for geodynamic processes is illustrated by two numerical experiments: i) volumetric-strain inducing phase changes in amagmatic slow-spreading ridges, and ii) subducting slabs.