  We consider a molecular dynamics method, the so-called flea gas for computing
the evolution of entanglement after inhomogeneous quantum quenches in an
integrable quantum system. In such systems the evolution of local observables
is described at large space-time scales by the Generalized Hydrodynamics
approach, which is based on the presence of stable, ballistically propagating
quasiparticles. Recently it was shown that the GHD approach can be joined with
the quasiparticle picture of entanglement evolution, providing results for
entanglement growth after inhomogeneous quenches. Here we apply the flea gas
simulation of GHD to obtain numerical results for entanglement growth. We
implement the flea gas dynamics for the gapped anisotropic Heisenberg XXZ spin
chain, considering quenches from globally homogeneous and piecewise homogeneous
initial states. While the flea gas method applied to the XXZ chain is not exact
even in the scaling limit (in contrast to the Lieb--Liniger model), it yields a
very good approximation of analytical results for entanglement growth in the
cases considered. Furthermore, we obtain the {\it full-time} dynamics of the
mutual information after quenches from inhomogeneous settings, for which no
analytical results are available.
