Variable-stiffness panels have previously shown enhanced buckling performance compared to constant-stiffness panels, due to the beneficial load and stiffness redistribution. Two semianalytical models based on the Galerkin method and the Ritz method have been developed in order to solve the buckling problem of variable-stiffness composite panels (plates and shallow cylindrical shells). The laminates considered are assumed to be symmetric, which results in zero extensionbending couplings. However, the bending-twisting couplings, D16 and D26, are retained in the model due to their dramatic influence on the convergence of the predicted buckling loads. In the first model the governing equations for composite plates and shallow cylindrical shells with variable stiffness are derived and then solved using the Galerkin method. The variable stiffness in this model is approximated using two-dimensional Fourier series, which, however, appears to be less accurate and less efficient as compared to the Ritz method. The Ritz method avoids using Fourier series to approximate the stiffness. Instead, the stiffness is exactly expressed in the integral of the energy functional. The detailed derivations of the energy functional which were rarely shown in literature have been presented. The buckling analysis of this model comprises two main steps. First, the in-plane loads are calculated by applying the principle of minimum complementary energy in pre-buckling state. Second, the critical buckling loads are determined from the stability equations which are obtained from the total energy functional through applying either adjacent equilibrium criterion or the principle of minimum potential energy. The total energy functional for stability analysis in this model is expressed in terms of out-of-plane displacement and the Airy stress function, which appears to be a combination of negative membrane complementary energy, bending strain energy and external work. In order to ensure fast convergence, several shape functions used in Ritz method were investigated. The in-plane loads were approximated either using the beam characteristic function or polynomial function that were orthogonalized by the Gram-Schmidt process. In addition, the predictions of in-plane loads using sine and cosine function that do not satisfy the boundary conditions of in-plane loads were significantly improved by using Lagrange multiplier method. The out-of-plane displacement was approximated using either sine function or orthogonalized polynomials. The influence of the shape functions on the convergence of the predicted buckling load was analysed and discussed for different examples. These examples included plates and shells with either constant-stiffness or variable stiffness, within which D16 and D26 were either zero or non-zero. The boundary conditions considered for these examples were four-edge simply-supported in current thesis. However, the developed model can be easily extended to consider other boundary conditions. The model can solve the buckling problem of variable-stiffness panels under prescribed inplane loads (¯Nx, ¯Ny, ¯Nxy) or prescribed in-plane displacements (¯u, ¯v). In the current thesis, only the prescribed loads ¯Nx, ¯Nxy and the prescribed end-shortenings ¯u have been investigated and compared to Abaqus model; all the results satisfactorily match the results of Abaqus models. In addition, the model have been proved to be able to predict the buckling loads of shallow cylindrical shells with variable curvatures.