We construct a Galerkin finite element method for the numerical approximation of weak solutions to a recent micro-macro bead-spring model for two-phase flow of dilute polymeric solutions derived by methods from nonequilibrium thermodynamics ([Grün, Metzger, M3AS 26 (2016) 823–866]). The model consists of Cahn-Hilliard type equations describing the evolution of the fluids and the unsteady incompressible Navier-Stokes equations in a bounded domain in two or three spatial dimensions for the velocity and the pressure of the fluids with an elastic extra-stress tensor on the right-hand side in the momentum equation which originates from the presence of dissolved polymer chains. The polymers are modeled by dumbbells subjected to a finitely extensible, nonlinear elastic (FENE) spring-force potential. Their density and orientation are described by a Fokker-Planck type parabolic equation with a center-of-mass diffusion term. We perform a rigorous passage to the limit as the spatial and temporal discretization parameters simultaneously tend to zero, and show that a subsequence of these finite element approximations converges towards a weak solution of the coupled Cahn-Hilliard-Navier-Stokes-Fokker-Planck system. To underline the practicality of the presented scheme, we provide simulations of oscillating dilute polymeric droplets and compare their oscillatory behaviour to the one of Newtonian droplets.