The Gibbs ensemble algorithm is implemented to determine the vapour-liquid phase coexistence of a purl fluid interacting via a two-body Lennard-Jones + three-body Axilrod-Teller potential. The contribution of both two-body and three-body interactions an calculated exactly. The results are compared with both experiment and two-body only simulation data, The position of the vapour branch of the coexistence curve is almost unaffected by the inclusion of three-body interactions. In contrast, the liquid branch occurs at substantially lower densities compared with Lennard-Jones simulation data. However, the approach to the critical point is improved by including three-body interactions, and the estimated critical point is in good agreement with the experiment.