The velocity distribution function, f1, for a (linear) shear flow of a system of rigid inelastically colliding disks in a plane is measured by applying a novel algorithm to results of (MD) simulations involving 200 000 particles. The need to consider such a relatively large system is explained. It is found that f1 is well fitted by an exponent of a second-order polynomial in the norm of the fluctuating velocities with angle-dependent coefficients (which also depend on the density and the granular temperature). Other characterizations of the system studied in this paper are presented as background material. A hitherto unnoticed property of systems with Lees–Edwards boundary conditions has been discovered and its origin is briefly explained.