The goal of the paper is to analyze a Gaudin model for a polynomial representation of the Kohno-Drinfeld Lie algebra associated with the multinomial distribution. The main result is the construction of an explicit basis of the space of polynomials consisting of common eigenfunctions of Gaudin operators in terms of Aomoto-Gelfand hypergeometric series. The construction shows that the polynomials in this basis are also common eigenfunctions of the operators for a dual Gaudin model acting on the degree indices, and therefore they provide a solution to a multivariate discrete bispectral problem.