An improved boundary element method for non-linear viscoelastic flow analysis is reported. In this method, the domain integral representing the non-linear effects is calculated using a more efficient approximating technique. This is achieved by first transforming the domain integral into a form that can be approximated by particular solutions to the original problem. These particular solutions are in fact expressed as a linear combination of radial basis functions, whose coefficients are found by data fitting technique. As a result, the numerical computation of the volume integral is eliminated and a significant reduction of CPU time is achieved. The routines are tested with simple flows and then applied to solve complex three-dimensional direct and inverse extrusion problems of polymeric fluids, such as thermoplastic melts. Inverse extrusion process, where an extrusion die profile needs to be computed for a given profile extrudate, is a very important practical engineering application which is successfully analysed by the present method. Relative to a previous BEM implementation where the volume integral is computed directly using numerical quadratures, a CPU time reduction ranging from 40 to 70% is achieved.