This paper is an extension to 3D FENE model of previous work [A. Lozinski, C. Chauviere, A fast solver for Fokker-Planck equation applied to viscoelastic flows calculations: 2D FENE model, J. Comput. Phys. 189 (2) (2003) 607-625] on the simulation of viscoelastic flows using a Fokker-Planck equation, where the 2D FENE model was considered. The Fokker-Planck equation here is discretized using spherical harmonics and Lagrange interpolating polynomial as basis functions. The scheme presented here exhibits the same attractive features as in [A. Lozinski, C. Chauviere, A fast solver for Fokker-Planck equation applied to viscoelastic flows calculations: 2D FENE model, J. Comput. Phys. 189 (2) (2003) 607-625], i.e. noise-free solutions and efficiency. The method is tested for the benchmark problem of a flow past a confined cylinder. Comparisons between the 2D FENE model and the 3D FENE model are made and it is found that the differences between those two models are small for moderate Deborah numbers. It is also found that stochastic simulations of 3D FENE model allow us to reach significantly higher Deborah numbers than that achievable by a similar Fokker-Planck simulations. (C) 2004 Elsevier B.V. All right reserved.