We propose a novel class of generalized skew-normal densities that improves the flexibility of empirical distributions and can systematically capture skewness, heavy tails, and multimodality. We extend the so-called flexible generalized skew-normal (FGSN) density developed by Y. Ma and M. G. Genton in 2004. The main novelty is the existence of a fifth-order degree term in the polynomial that appears in the cumulative distribution function of such a density. In this case, we prove that our density has at most three modes under certain conditions for the parameters. Leveraging this new approach eases the modeling of data consisting of three subpopulations. For validation, we present examples of both univariate and multivariate cases using demographic data from the Human Mortality Data Base (HMD).