Simple analytical and finite element models are widely employed by practising engineers for the stress analysis of beam structures, because of their simplicity and acceptable levels of accuracy. However, the validity of these models is limited by assumptions of material heterogeneity, geometric dimensions and slenderness, and by Saint-Venant’s Principle, i.e. they are only applicable to regions remote from boundary constraints, discontinuities and points of load application. To predict accurate stress fields in these locations, computationally expensive three-dimensional (3D) finite element analyses are routinely performed. Alternatively, displacement based high-order beam models are often employed to capture localised three-dimensional stress fields analytically. Herein, a novel approach for the analysis of beam-like structures is presented. The approach is based on the Unified Formulation by Carrera and co-workers, and is able to recover complex, 3D stress fields in a computationally efficient manner. As a novelty, purposely adapted, hierarchical polynomials are used to define cross-sectional displacements. Due to the nature of their properties with respect to computational nodes, these functions are known as Serendipity Lagrange polynomials. This new cross-sectional expansion model is benchmarked against traditional finite elements and other implementations of the Unified Formulation by means of static analyses of beams with different complex cross-sections. It is shown that Serendipity Lagrange elements solve some of the shortcomings of the most commonly used Unified Formulation beam models based on Taylor and Lagrange expansion functions. Furthermore, significant computational efficiency gains over 3D finite elements are achieved for similar levels of accuracy.