First BGG operators are a large class of overdetermined linear differential operators intrinsically associated to a parabolic geometry on a manifold. The corresponding equations include those controlling infinitesimal automorphisms, higher symmetries and many other widely studied PDE of geometric origin. The machinery of BGG sequences also singles out a subclass of solutions called normal solutions. These correspond to parallel tractor fields and hence to (certain) holonomy reductions of the canonical normal Cartan connection. Using the normal Cartan connection, we define a special class of local frames for any natural vector bundle associated to a parabolic geometry. We then prove that the coefficient functions of any normal solution of a first BGG operator with respect to such a frame are polynomials in the normal coordinates of the parabolic geometry. A bound on the degree of these polynomials in terms of representation theory data is derived. For geometries locally isomorphic to the homogeneous model of the geometry we explicitly compute the local frames mentioned above. Together with the fact that on such structures all solutions are normal, we obtain a complete description of all first BGG solutions in this case. Finally, we prove that in the general case the polynomial system coming from a normal solution is the pull-back of a polynomial system that solves the corresponding problem on the homogeneous model. Thus we can derive a complete list of potential normal solutions on curved geometries. Moreover, questions concerning the zero locus of solutions, as well as related finer geometric and smooth data, are reduced to a study of polynomial systems and real algebraic sets.