The problem of development of the high-order numerical scheme for 2D flows simulation in fully Lagrangian meshless vortex methods is considered. The main computation variable is vorticity; its evolution in the flow (transfer and diffusion) is simulated through the motion of vortex elements taking into account the so-called diffusive velocity according to the Viscous Vortex Domains method. Vorticity generation on the airfoil surface line is described by the integral equation, for which solution the original numerical scheme is developed. It deals with piecewise-quadratic vorticity distribution along the airfoil surface line and provides the third order of accuracy. However, in the framework of traditional approach to solution of boundary integral equation, in order to take into account the influence of nearby vortex elements, it is necessary to refine surface line discretization. In practice, it is hardly possible to achieve surface mesh refinement, sufficient for correct computation of the contributions of the vortex elements in the flow, so semi-analytical correction procedure is suggested. All computational formulae are derived.