The form and evolution of multi-phase biomembranes is of fundamental impor-tance in order to understand living systems. In order to describe these membranes,we consider a mathematical model based on a Canham–Helfrich–Evans two-phaseelastic energy, which will lead to fourth order geometric evolution problems involv-ing highly nonlinear boundary conditions. We develop a parametric finite elementmethod in an axisymmetric setting. Using a variational approach, it is possible toderive weak formulations for the highly nonlinear boundaryvalue problems suchthat energy decay laws, as well as conservation properties,hold for spatially discre-tised problems. We will prove these properties and show thatthe fully discretisedschemes are well-posed. Finally, several numerical computations demonstrate thatthe numerical method can be used to compute complex, experimentally observedtwo-phase biomembranes.