A finite element method is presented for a nonlocal Timoshenko beam model recently proposed by the authors. The model relies on the key idea that nonlocal effects consist of long-range volume forces and moments exchanged by non-adjacent beam segments, which contribute to the equilibrium of a beam segment along with the classical local stress resultants. The long-range volume forces/moments are linearly depending on the product of the volumes of the interacting beam segments, and their relative motion measured in terms of the pure beam deformation modes, through appropriate attenuation functions governing the spatial decay of nonlocal effects. In this paper, the beam model is reformulated within a variational framework involving a consistent total elastic potential energy functional. The latter serves as a basis to derive a suitable finite element formulation of the equilibrium equations. A local stiffness matrix and a nonlocal stiffness matrix contribute to the global stiffness matrix. While the local stiffness matrix is obtained by a standard assemblage of the classical element stiffness matrices, the nonlocal stiffness matrix is built as the sum of component matrices, each involving the stiffness of the long-range interactions between a couple of finite elements. A remarkable result is that, for most common attenuation functions of nonlocal effects, exact closed-form solutions can be found for every element of the nonlocal stiffness matrix. Numerical applications are presented for a variety of nonlocal parameters, including a comparison with experimental data.