This talk focuses on improving the robustness of the functional subspace built using Proper Orthogonal Decomposition (POD). Since a POD basis is able to give an optimal representation of the kinetic energy included in the snapshots database generated with some given input parameters, this same basis is not adapted to represent flow dynamics generated with other input parameters. Our aim is thus to build a robust POD basis that accurately represents the solution over a desired input parameter subspace by enlarging the statistical database. We present a systematic method to sample the input parameter subspace. The basic idea is to add to the existing database, snapshots of the solution for which the POD approximation error is maximal. Then the new POD basis is obtained from the set of snapshots corresponding to such different dynamics. This is the so called Greedy sampling. The approach we follow is similar: it is based on finding the centroid of a region around the point where an estimate of the POD approximation error is maximal. We show a priori numerical evidence that the Navier-Stokes residuals are a reliable estimate of the POD approximation error. Results relative to a 2D confined square cylinder wake flow are presented. The input parameter subspace is represented by an interval of Reynolds numbers that corresponds to periodical laminar flows. We show that a judicious choice of the sampling Reynolds numbers leads to a POD basis that minimizes the average approximation error on the chosen interval.