In this article, two analytical solutions, one semi-analytical solution, and one numerical solution are proposed to calculate the global critical buckling load of buildings employing a continuum model parallel coupling two Timoshenko beams. The parallel coupling of two Timoshenko beams represents the four types of building behavior: global bending, global shear, local bending, and local shear, the latter of which has been widely ignored in the literature. The model is associated with one translational kinematic field and two rotational fields. The equilibrium equations, constitutive laws, and boundary conditions are derived using a variational approach by applying Hamilton’s principle to the relevant Lagrangian function. Specifically, to solve the classical problem of buildings with uniform properties along their height, the first analytical solution proposes exact closed-form equations for a vertical load applied at the top and a closed-form analytical solution based on the decomposition of two subsystems for a polynomial vertical load uniformly applied along the height of the beam. Meanwhile, the semi-analytical solution exactly solves the differential equation associated with a polynomial vertical load profile and is based on coefficient iteration, converging with only two or three iterations. To address the challenge of buildings with variable properties along their height, the combined application of the continuous method and the transfer matrix method allows us to propose an efficient numerical method leading to a 6 × 6 transfer matrix, constant regardless of the number of discretizations and derived analytically, which results in a computational advantage by drastically reducing computational costs. The numerical results are encouraging and show good agreement with the results of the finite element method, suggesting that the proposed method can be used for engineering purposes in both the preliminary and final phases of the structural analysis of tall buildings.