We propose a method of constructing analytical, closed-form expressions for electrostatic/Newtonian potentials of non-uniform polyhedral bodies, in which the density distributions are polynomials of coordinates. Possible applications of the proposed method are spread from astronomy to nanotechnology. The method is based on the use of the generating function for the potential. Explicit expressions for the potential are derived in the case of quadratic or cubic coordinate dependence of the density within a polyhedral body.