This talk is devoted to the numerical approximation of the bidimensional bitemperature Euler system. This model is a nonconservative hyperbolic system describing an out of equilibrium plasma in a quasi-neutral regime. This system is non conservative because it involves products between velocity and pressure gradients that cannot be transformed into a divergential form. We develop a second order numerical scheme by using a discrete BGK relaxation model. The second order extension is based on a subdivision of each cartesian cell into four triangles to perform affine reconstructions of the solution. Such ideas have been developed in the litterature for systems of conservation laws. We show here how they can be used in our nonconservative setting. The numerical method is implemented and tested in the last part of the paper.