A backward differences formulae (BDF) scheme, is proposed to simulate the deformation of a viscous incompressible Newtonian fluid domain in time, which is driven solely by the boundary curvature. The boundary velocity field of the fluid domain is obtained by writing the governing Stokes equations in terms of an integral formulation that is solved by a boundary element method. The motion of the boundary is modelled by considering the boundary curve as material points. The trajectories of those points are followed by applying the Lagrangian representation for the velocity. Substituting this representation into the discretized version of the integral equation yields a system of non-linear ODEs. Here the numerical integration of this system of ODEs is outlined. It is shown that, depending on the geometrical shape, the system can be stiff. Hence, a BDF-scheme is applied to solve those equations. Some important features with respect to the numerical implementation of this method are high-lighted, like the approximation of the Jacobian matrix and the continuation of integration after a mesh redistribution. The usefulness of the method for both two-dimensional and axisymmetric problems is demonstrated.