In this paper, based on the nonlocal Timoshenko beam theory, a nonlinear model is presented for the vibrational behavior of carbon nanotubes (CNTs) embedded in elastic medium in thermal environment. Using the Timoshenko beam theory and nonlocal elasticity of Eringen, the influences of rotary inertia, transverse shear deformation and small scale effect are taken into account. To model the interaction forces between walls, whether adjacent or non-adjacent, the van der Waals interlayer interactions are considered. The harmonic balance method (HBM) is used for the solution of the set of nonlinear governing equations and the frequency function of the system for the simply-supported boundary conditions is derived. Compared to the incremental harmonic balance method which has been employed in the previous studies, the HBM is simpler and has a reasonable accuracy. The effects of geometrical parameters of nanotubes such as the number of walls, the ratio of length to outer diameter and environmental conditions such as elastic medium modulus, temperature and also the effect of nonlocal parameter on the nonlinear frequency are investigated. The presented nonlinear vibration analysis is of a general form, so that they are applicable for CNTs with arbitrary number of walls. The obtained results for single-, double- and triple-walled CNTs indicate that with an increase in the number of walls, elastic medium modulus, aspect ratio and temperature, the value of nonlinear frequency tends to that of its linear counterpart. Also, a comparison between the results of the Timoshenko beam theory and those of Euler-Bernoulli beam theory shows that the difference between the frequency responses of these theories is significant for short CNTs, but, as the length increases, the difference between the results becomes negligible.