  We consider systematic numerical approximation of a viscoelastic phase
separation model that describes the demixing of a polymer solvent mixture. An
unconditionally stable discretisation method is proposed based on a finite
element approximation in space and a variational time discretization strategy.
The proposed method preserves the energy-dissipation structure of the
underlying system exactly and allows to establish a fully discrete nonlinear
stability estimate in natural norms based on the concept of relative energy.
These estimates are used to derive order optimal error estimates for the method
under minimal smoothness assumptions on the problem data, despite the presence
of various strong nonlinearities in the equations. The theoretical results and
main properties of the method are illustrated by numerical simulations which
also demonstrate the capability to reproduce the relevant physical effects
observed in experiments.
