•Three-dimensional solution of functionally graded shells is analyzed.•Differential Quadrature Method (DQM) is used.•Chebyshev-Gauss-Lobatto grid and Lagrange interpolation polynomials are used.•Shells are subjected to bi-sinusoidal and uniform distributed load.•Benchmark results are provided. A numerical solution for the three-dimensional static analysis of functionally graded shells with constant curvature is presented. The solution is based on three-elasticity equations written in orthogonal curvilinear coordinates which are valid for spherical, cylindrical shell panels and rectangular plates. The equations in term of the mid-surface variables are solved using a summation of harmonics in term of Navier method which is valid only for simply supported structures. The equations in term of the thickness direction are solved numerically by the Differential Quadrature method (DQM) which permitted to easily calculate the approximate derivative of a function using a weighting sum of the functions evaluated in a certain grid. The layers of the structure are discretized separately by the Chebyshev-Gauss-Lobatto grid and Lagrange interpolation polynomials are considered as the basis functions. The inter-laminar continuity of transverse shear is imposed as part of the boundary conditions for the presented method. The boundary conditions of out-of-plane stresses at the top and the bottom due to the applied loads on the shell are also considered for the analysis, as a result this method can predict the correct behavior of through-the-thickness distribution of transverse stresses. This method permitted easily to discretize the material in term of the thickness direction and several types of single functionally graded layer and sandwich structures with functionally graded core are analyzed. Several shells subjected to bisinusoidal and uniform distributed load are analyzed. The results are compared with other three dimensional solutions proposed in the literature and accurate two dimensional models.