  The upcoming generation of galaxy surveys will probe the distribution of
matter in the universe with unprecedented accuracy. Measurements of the matter
power spectrum at different scales and redshifts will provide stringent
constraints on the cosmological parameters. However, on non-linear scales this
will require an accurate evaluation of the covariance matrix. Here, we compute
the covariance matrix of the 3D matter density power spectrum for the
concordance $\Lambda$CDM cosmology from an ensemble of N-body simulations of
the Dark Energy Universe Simulation - Parallel Universe Runs (DEUS-PUR). This
consists of 12288 realisations of a $(656\,h^{-1}\,\textrm{Mpc})^3$ simulation
box with $256^3$ particles. We combine this set with an auxiliary sample of 96
simulations of the same volume with $1024^3$ particles. We find N-body mass
resolution effect to be an important source of systematic errors on the
covariance at high redshift and small intermediate scales. We correct for this
effect by introducing an empirical statistical method which provide an accurate
determination of the covariance matrix over a wide range of scales including
the Baryon Oscillations interval. Contrary to previous studies that used
smaller N-body ensembles, we find the power spectrum distribution to
significantly deviate from expectations of a Gaussian random density field at
$k\gtrsim 0.25\,h\,\textrm{Mpc}^{-1}$ and $z<0.5$. This suggests that in the
case of finite volume surveys an unbiased estimate of the ensemble averaged
band power at these scales and redshifts may require a careful assessment of
non-Gaussian errors more than previously considered.
