An integrated nonlinear couple stress-surface energy continuum model is developed to study the nonlinear vibration characteristics of size-dependent functionally graded nanobeams for the first time. The nanobeam theory is formulated based on the Timoshenko kinematics, augmented by von Kármán’s geometric nonlinearity. The modified couple stress and Gurtin–Murdoch surface elasticity theories are incorporated to capture the long-range interaction and surface energy, respectively. Unlike existing Timoshenko nanobeam models, the effects of surface elasticity, residual surface stress, surface mass density and Poisson’s ratio, in addition to bending and axial deformations, are incorporated in the newly developed model. A power law function is used to model the material distribution through the thickness of the beam, considering the gradation of bulk and surface material parameters. A variational formulation of the nonlinear nonclassical governing equations and associated nonclassical boundary conditions is established by employing Hamilton’s principle. The generalized differential quadrature method is exploited in conjunction with either the Pseudo-arclength continuation or Runge–Kutta method to solve the problem with an exact implementation of the nonclassical boundary conditions. The formulation and solution procedure presented are verified by comparing the obtained results with available ones. Based on the parametric study, it is concluded that the nonclassical boundary conditions, material length scale parameter, residual surface stress, surface elasticity, bulk elasticity modulus, gradient index, nonlinear amplitude and thickness have important influences on the linear and nonlinear vibration responses of functionally graded Timoshenko nanobeams.