A new formulation of the initial-boundary value problem of the N-th order orthotropic shells theory based on the methods of analytical mechanics of continuum systems is proposed below. Some smooth base surface is introduced, and curvilinear coordinates are determined on the two-dimensional manifold corresponding to the base surface. The model of a non-thin shell as a three-dimensional elastic body is defined by a set of field variables of the first kind, and the spatial and boundary densities of the Lagrange functional. As a field variable of the first kind, a pseudovector of displacement is considered, given by covariant components of the vector of true displacement in the basis of the tangent bundle of a two-dimensional manifold. The components of the generalized stress pseudovector on sections with normal unit co-directed to the base vector of one of the curvilinear coordinates are defined by differentiating the spatial density of the Lagrange functional by covariant derivatives of the components of the displacement pseudovector in the selected direction. By Legendre transformation of the spatial density of the Lagrange functional with respect to to the covariant derivatives of the field variables, the spatial density of the mixed functional is obtained, depending on the state variables – the introduced stress pseudovector, displacement and velocity pseudovectors, as well as covariant derivatives of the components of the displacement pseudovector in the second coordinate direction. A biorthogonal base system of functions of the dimensionless normal coordinate is introduced, a system of new state variables defined on the tangent bundle of a two-dimensional manifold corresponding to the basis surface is determined by the series coefficients of state variables with respect to the biorthogonal base system. The surface and contour densities of the mixed functional defining the two-dimensional shell model are obtained by the corresponding reduction of the spatial and boundary densities of the mixed functional while retaining first N+1 series coefficient. The Euler-Lagrange equations, which follow from the Hamilton-Ostrogradsky principle, are resolved hence with respect to covariant derivatives of new state variables and, in a certain sense, are equivalent to the known Routh equations of the analytical mechanics of discrete systems. The application of the generalized Routh equations of the N-th order shell theory to problems of normal wave dispersion leads to a spectral problem linear with respect to the wavenumber, which allows one to construct as well real as imaginary and complex branches of dispersion curves corresponding to propagating and evanescent wave modes. A solution of the spectral problem for an isotropic plane layer is obtained and the solution convergence based on the N-th order theory to the exact solution in the second quadrant of the complex plane at some nodes of the Mindlin grid is studied for imaginary branches.