A powerful extrapolation scheme is proposed to determine the vapour-liquid equilibrium (VLE) curves of pure molecular fluids over a finite temperature range by performing a single Gibbs ensemble Monte Carlo (GEMC) simulation. The VLE curves of the various thermodynamic quantities are extrapolated as functions of the temperature by second-order Taylor series. The coefficients of the Taylor series, which are the temperature derivatives of these quantities along the VLE curves, can be calculated from a single GEMC simulation using Clapeyron-like equations and fluctuation formulas containing double or triple correlations. It is shown that the application of a Padé approximant considerably widens the temperature range where the extrapolation is accurate. The efficiency of the method is demonstrated by applying it on the Lennard-Jones fluid and a three-site model of carbon dioxide. The procedure is found to be a very useful tool that is able to reproduce the VLE curves over a considerable temperature range. The procedure makes it possible to calculate the saturation heat capacity from fluctuation formulas.