A (3,2) high order zigzag beam element based on unified beam theory is developed. A unified zigzag function family is derived based on the global displacement and the initial rigidity of the layup. In the determination of zigzag functions, two approaches are proposed. The first is to solve a linear system, where the matrix’s rank is 1 order less than the size. The second is to derive a recursive solution to determine the zigzag functions. Any specific beam theories could be embedded here. The traction free condition on top/bottom surfaces is not used in this paper, instead, a more general traction condition on top/bottom surfaces is utilized. In the axial direction, the Consistent Orthogonal Basis Function Space is applied, such that the basis functions are very identical to mode shape functions. The finite element equation implementations are also presented.In the numerical tests, several bench problems are solved to verify the accuracy of the method. It is observed that the method has a fast convergence rate for displacement and its derivatives. Then, the laminated beam with sharp change of Young’s modulus along thickness direction is also tested. Finally, the shear stress of beam under tangential traction is studied.