We devise a fully self-consistent simulation pipeline for the first time to study the interaction between dark matter and dark energy. We perform convergence tests and show that our code is accurate on different scales. Using the parameters constrained by Planck, Type Ia Supernovae, Baryon Acoustic Oscillations (BAO) and Hubble constant observations, we perform cosmological N-body simulations. We calculate the resulting matter power spectra and halo mass functions for four different interacting dark energy models. In addition to the dark matter density distribution, we also show the inhomogeneous density distribution of dark energy. With this new simulation pipeline, we can further refine and constrain interacting dark energy models.