Thick composite-layered shells under hydrostatic pressure are analyzed for buckling and postbuckling response using axisymmetric solid elements. The numerical approach is based on p-version finite elements in conjunction with appropriate trigonometric functions. Particular attention is given to the evaluation of interlaminar stresses in the postbuckling range. The postbuckling response is determined using an asymptotic approach. It is found that layered shells with Z (the modified Batdorf parameter) < 1000 are imperfection-sensitive with respect to local buckling; shells with relatively small radius to thickness ratios develop significant interlaminar shearing stresses, of the order of 500% of the buckling pressure, for buckling displacements of the order of one-tenth of the shell thickness. Together, the imperfection-sensitivity and the high interlaminar stresses can precipitate failure before the classical critical pressures are reached. These factors must be considered along with the potential for interaction of local and overall buckling in optimally designed, stiffened shell configurations. Numerical analyses demonstrate the effectiveness of the p-version approach in achieving rapid convergence of the buckling pressures, the imperfection-sensitivity parameter and the magnitude of interlaminar stresses. Wherever possible, comparisons have been made with results available currently in the literature. It is found that for long shells (Z > 1000), the Donnell type formulations can seriously underestimate the buckling pressures. The present formulation based on three-dimensional nonlinear elasticity offers an accurate, powerful and yet computationally effective means of analysis of thick shells.