This work presents the analytical and numerical formulation for the analysis of thick spherical shells, under internal and external pressure loads, utilizing the Finite Element Method with unidimensional elements, containing four nodal points each. The numerical model incorporates the simplifications resulting from the axial symmetry, which determined the choice of the shape of the elements. Some boundary conditions of the elements allow the successful representation of the geometry of the chosen shell. The first refers to the natural boundary condition of zero shear stress at the internal and external surfaces of the shell. The second condition aims to guarantee the continuity between two adjacent elements. Finally, the third condition will represent the chosen fixation conditions of an element. The first condition is guaranteed through the imposition of the linking between the model’s degrees-of-freedom, while the other two will be guaranteed numerically by using the Penalty Method. Displacements of the model’s nodal points are represented by two displacement fields, longitudinal and radial. A total of six degrees-of-freedom per nodal point are utilized for representing these displacement fields. The numerical model was implemented, and the obtained results were compared with existing analytical solutions, in order to ascertain the validity of the adopted methodology.