A high-order theory is developed for the analysis of beams with general mono-symmetric cross-sections. The theory represents the nonlinear distribution of the longitudinal normal stress across the section depth by a polynomial series expansion up to any order as specified by the analyst. The corresponding shear and transverse normal stresses are obtained by satisfying the 2D infinitesimal stress-flow equilibrium conditions. The resulting statically admissible stress fields are then used in conjunction with the principle of stationary complementary strain energy to formulate the governing compatibility equations and boundary conditions. Closed form solutions are then developed for general loading and boundary conditions. Comparisons with the theory of elasticity and 3D finite element analysis predictions showcase the ability of the present theory to naturally capture shear deformation effects, transverse normal stress effects, nonlinear longitudinal normal stress distributions in deep beams, as well as the effect of support height. Unlike conventional beam solutions that are based on postulated kinematic assumptions, which tend to converge to the displacement response from below, the present theory avoids introducing any kinematic assumptions and is shown to converge to the solution from above. The theory is applicable to beams with doubly symmetric or mono-symmetric cross-sections, with isotropic or orthotropic materials, and subjected to general loading and boundary conditions. The theory is shown to offer advantages compared to other theories when modelling deep beams and/or beams with supports that are offset from the centroidal axis.