We study variational inequality problems (VIPs) with involved mappings and feasible sets characterized by polynomial functions (namely, polynomial VIPs). We propose a numerical algorithm for computing solutions to polynomial VIPs based on Lagrange multiplier expressions and the Moment-SOS hierarchy of semidefinite relaxations. We also extend our approach to finding more or even all solutions to polynomial VIPs. We show that the method proposed in this paper can find solutions or detect the nonexistence of solutions within finitely many steps, under some general assumptions. In addition, we show that if the VIP is given by generic polynomials, then it has finitely many Karush-Kuhn-Tucker points, and our method can solve it within finitely many steps. Numerical experiments are conducted to illustrate the efficiency of the proposed methods.