We construct a general family of Galerkin methods for the numerical approximation of weak solutions to a coupled microscopic-macroscopic bead-spring model that arises from the kinetic theory of dilute solutions of polymeric liquids with noninteracting polymer chains. The model consists of the unsteady incompressible Navier-Stokes equations in a bounded domain Ω in R d, d=2 or 3, for the velocity and the pressure of the fluid, with an elastic extra-stress tensor as right-hand side in the momentum equation. The extra-stress tensor stems from the random movement of the polymer chains and is defined through the associated probability density function which satisfies a Fokker-Planck type parabolic equation, a crucial feature of which is the presence of a centre-of-mass diffusion term. We focus on finitely-extensible nonlinear elastic, FENE-type, dumbbell models. In the case of a corotational drag term we perform a rigorous passage to the limit as the spatial and temporal discretization parameters tend to zero, and show that a (sub)sequence of numerical solutions converges to a weak solution of this coupled Navier-Stokes-Fokker-Planck system.