In this talk, we present a strategy to numerically solve, via structure preserving high order Finite Volume (FV) and Discontinuous Galerkin (DG) methods, i) the general relativistic hydrodynamics equations (Euler), ii) the first order hyperbolic Z4 re-formulation of the Einstein field equations of 3+1 general relativity, as well as iii) the fully coupled Einstein + Euler system. Indeed, thanks to our hyperbolic reformulation of the Eistein field equations, we can study the dynamics of the matter coupled with the dynamics of the space-time, described by the theory of general relativity, under a unique formalism, without the need of using different numerical methods for matter and metric. From the numerical point of view, we make use of well-balanced techniques to enhance the resolution close to equilibrium profiles, a covariant GLM cleaning technique to assure the preservation of nonlinear involutions, the a posteriori FV limiter to avoid oscillations in the DG methods close to shock discontinuities, and a new filter for the conversion from the conserved to the primitive variables preventing the emergence of superluminal velocities near vacuum. We will close the talk by showing some challenging numerical results we have obtained, as the stable long term simulations of stationary black holes, including Kerr black holes with extreme spin, the evolution of a TOV star under perturbation in pure vacuum, and the head on collision of two punctures black holes.