  We present a new, monolithic first--order (both in time and space) BSSNOK
formulation of the coupled Einstein--Euler equations. The entire system of
hyperbolic PDEs is solved in a completely unified manner via one single
numerical scheme applied to both the conservative sector of the matter part and
to the first--order strictly non--conservative sector of the spacetime
evolution. The coupling between matter and space-time is achieved via algebraic
source terms. The numerical scheme used for the solution of the new monolithic
first order formulation is a path-conservative central WENO (CWENO) finite
difference scheme, with suitable insertions to account for the presence of the
non--conservative terms. By solving several crucial tests of numerical general
relativity, including a stable neutron star, Riemann problems in relativistic
matter with shock waves and the stable long-time evolution of single and binary
puncture black holes up and beyond the binary merger, we show that our new
CWENO scheme, introduced two decades ago for the compressible Euler equations
of gas dynamics, can be successfully applied also to numerical general
relativity, solving all equations at the same time with one single numerical
method. In the future the new monolithic approach proposed in this paper may
become an attractive alternative to traditional methods that couple central
finite difference schemes with Kreiss-Oliger dissipation for the space-time
part with totally different TVD schemes for the matter evolution and which are
currently the state of the art in the field.
