We study a generalization of Eulerian polynomials to the multivariate setting introduced by Brändén. Although initially these polynomials were introduced using the language of hyperbolic and stable polynomials, we manage to translate some restrictions of these polynomials to our real zero setting. Once we are in this setting, we focus our attention on the rigidly convex sets (RCSs) defined by these polynomials. In particular, we study the corresponding rigidly convex sets looking at spectrahedral relaxations constructed through the use of monic symmetric linear matrix polynomials (MSLMPs) of small size and depending polynomially (actually just cubically) on the coefficients of the corresponding polynomials. We analyze how good are the obtained spectrahedral approximations to these rigidly convex sets. We do this analysis by measuring the behavior along the diagonal, where we precisely recover the original univariate Eulerian polynomials. Thus we conclude that, measuring through the diagonal, our relaxation-based spectrahedral method for approximation of the rigidly convex sets defined by multivariate Eulerian polynomials is highly accurate. In particular, we see that this relaxation-based spectrahedral method for approximation of the rigidly convex sets defined by multivariate Eulerian polynomials provides bounds for the extreme roots of the corresponding univariate Eulerian polynomials that are better than these already found in the literature. All in all, this tells us that, at least close to the diagonal, the global outer approximation to the rigidly convex sets provided by this relaxation-based spectrahedral method is itself highly accurate.