We revisit in one dimension the waterbag method to solve numerically Vlasov-Poisson equations. In this approach, the phase-space distribution function $f(x,v)$ is initially sampled by an ensemble of patches, the waterbags, where $f$ is assumed to be
constant. As a consequence of Liouville theorem it is only needed to follow the evolution of the border of these waterbags, which can be done by employing an orientated, self-adaptive polygon tracing isocontours of $f$. This method, which is entropy conserving in essence, is very accurate and can trace very well non linear instabilities as illustrated by specific examples. As an application of the method, we generate an ensemble of single waterbag simulations with decreasing thickness, to perform a convergence study to the cold case. Our measurements show that the system relaxes to a steady state where the gravitational potential profile is a power-law of slowly varying index $beta$, with $beta$ close to $3/2$ as found in the literature. However, detailed analysis of the properties of the gravitational potential shows that at the center, $beta > 1.54$. Moreover, our measurements are consistent with the value $beta=8/5=1.6$ that can be analytically derived by assuming that the average of the phase-space density per energy level obtained at crossing times is conserved during the mixing phase. These results are incompatible with the logarithmic slope of the projected density profile $beta-2 simeq -0.47$ obtained recently by Schulz et al. (2013) using a $N$-body technique. This sheds again strong doubts on the capability of $N$-body techniques to converge to the correct steady state expected in the continuous limit.
