A three-dimensional boundary-integral algorithm is used to study thermocapillary interactions of two deformable drops in the presence of bulk-insoluble, non-ionic surfactant. The primary effect of deformation is to slow down the rate of film drainage between drops in close approach and prevent coalescence in the absence of van der Waals forces. Both linear and non-linear models are used to describe the relationship between interfacial tension and surfactant surface concentration. In the linear model, non-monotonic behavior of the minimum separation between the drops as a function of the surface Peclet number Pe_s is observed for equal drop and external medium viscosities and thermal conductivities. For bubbles with zero drop-to-medium viscosity and thermal conductivity ratios, however, the minimum separation increases with Pe_s. There is a nearly linear relationship between the minimum drop separation and elasticity E. In the simplest non-linear equation of state, the product of the temperature and the surfactant concentration is retained by allowing non-zero values of the dimensionless gas constant Λ. For Λ=O(0.05), it is possible for the smaller drop to move faster than the larger drop. In the Langmuir adsorption framework, the tendency of the smaller drop to catch up to the larger one decreases as the ratio of the equilibrium to maximum surfactant surface concentration increases. Finally, in the Frumkin model, a minimum in the drop separation occurs as a function of the interaction parameter λ_F for trajectories with all other parameters held constant.