Numerical method is presented for simulation of the deformation, drainage and rupture of axisymmetric film (gap) between colliding drops in the presence of film soluble surfactants under the influence of van der Waals forces at small capillary and Reynolds numbers and small surfactant concentrations. The mathematical model is based on the lubrication equations in the gap between drops and the creeping flow approximation of Navier-Stokes equations in the drops, coupled with velocity and stress boundary conditions at the interfaces. A non-uniform surfactant concentration on the interfaces, related with that in the film, leads to a gradient of the interfacial tension which in turn leads to additional tangential stress on the interfaces (Marangoni effects). Both film and interface surfactant concentrations, related via adsorption isotherm, are governed by a convection-diffusion equation. The numerical method consists of: Boundary integral method for the flow in the drops; Finite difference method for the flow in the gap, the position of the interfaces and the surfactant concentration on the interfaces, as well as in the film. Second order approximation of the spatial terms on adaptive non-uniform mesh is constructed in combination with Euler explicit scheme for the time discretization. For the convection-diffusion equation in the film first order implicit and Crank-Nicolson time integration schemes are used as well. Tests and comparisons are performed to show the accuracy and stability of the presented numerical method.