The shapes of vesicles formed by lipid bilayers with phase separation are governed by a bending energy with phase dependent material parameters together with a line energy associated with the phase interfaces. We present a numerical method to approximate solutions to the Euler-Lagrange equations featuring triangulated surfaces, isoparametric quadratic surface finite elements and the phase field approach for the phase separation. Furthermore, the method involves an iterative solution scheme that is based on a relaxation dynamics coupling a geometric evolution equation for the membrane surface with a surface Allen-Cahn equation. Remeshing and grid adaptivity are discussed, and in various simulations the influence of several physical parameters is investigated.