We consider the problem of refining a parameter set to ensure that the behaviors of a dynamical system satisfy a given property. The dynamics are defined through parametric polynomial difference equations and their Bernstein representations are exploited to enclose reachable sets into parallelotopes. This allows us to achieve more accurate reachable set approximations with respect to previous works based on axis-aligned boxes. Moreover, we introduce a symbolical precomputation that leads to a significant improvement on time performances. Finally, we apply our framework to some epidemic models verifying the strength of the proposed method.