A nonlocal Bernoulli–Euler p-version finite-element (p-FE) is developed to investigate nonlinear modes of vibration and to analyze internal resonances of beams with dimensions of a few nanometers. The time domain equations of motion are transformed to the frequency domain via the harmonic balance method (HBM), and then, the equations of motion are solved by an arc-length continuation method. After comparisons with published data on beams with rectangular cross section and on carbon nanotubes (CNTs), the study focuses on the nonlinear modes of vibration of CNTs. It is verified that the p-FE proposed, which keeps the advantageous flexibility of the FEM, leads to accurate discretizations with a small number of degrees-of-freedom. The first three nonlinear modes of vibration are studied and it is found that higher order modes are more influenced by nonlocal effects than the first mode. Several harmonics are considered in the harmonic balance procedure, allowing us to discover modal interactions due to internal resonances. It is shown that the nonlocal effects alter the characteristics of the internal resonances. Furthermore, it is demonstrated that, due to the internal resonances, the nonlocal effects are still noticeable at lengths that are longer than what has been previously found.