This work is dedicated to the construction of numerical techniques for the models of viscoelastic fluids that result from polymer kinetic theory. Our main contributions are as follows: Inspired by the interpretation of the Oldroyd B model of dilute polymer solutions as a suspension of Hookean dumbbells in a Newtonian solvent, we have constructed new numerical methods for this model that respect some important properties of the underlying differential equations, namely the positive definiteness of the conformation tensor and an energy estimate. These methods have been implemented on the basis of a spectral discretization for simple Couette and Poiseuille planar flows as well as flow past a cylinder in a channel. Numerical experiments confirm the enhanced stability of our approach. Spectral methods have been designed and implemented for the simulation of mesoscopic models of polymeric liquids that do not possess closed-form constitutive equations. The methods are based on the Fokker-Planck equations rather than on the equivalent stochastic differential equations. We have considered the FENE dumbbell model of dilute polymer solutions and the Öttinger reptation model of concentrated polymer solutions. The comparison with stochastic simulation techniques has been performed in the cases of both homogeneous flows and the flow past a cylinder in a channel. Our method turned out to be more efficient in most cases.