In a recent study a geometrically nonlinear two-dimensional formulation of a Euler-Bernoulli beam element was presented. It uses an integrated form of equilibrium equations that are combined with the kinematic and generalized material equations, resulting in a set of three first-order differential equations. Finite differences are used to discretize the equations and transform the boundary value problem into an initial value problem using a technique inspired by the shooting method. In this work we extend the formulation to include transverse shear strain with the aim of modelling metamaterials with internal microstructure that corresponds to a two-dimensional lattice, manufactured using short beams or where shear strain variations through the thickness can be significant. The study focuses on stability analysis for critical load prediction. In particular, the implementation includes both the Reissner/Haringx and Engesser model. The validation is carried out through comparison with standard beam elements formulated in the finite-strain framework and with analytical solutions. This extension allows an increase in the range of choice of the specific metamaterial with user-defined properties. Moreover, the overall results achieved with the proposed computational approach can be used as a reference for the design requirements, enabling the description of a broad range of mechanical phenomena.