The aim of our work is the numerical modeling of two dimensional mechanical–thermal material mixing observed in stir welding process using a high order algorithm. This algorithm is based on coupling a meshless approach, a time discretization, a homotopy transformation, a development in Taylor series and a continuation method. The performance of the proposed model is the consideration of large deformations in the formulation of the posed problem. For the spatial treatment, we use the moving least squares approximation which will be applied directly to the strong form formulation of conservation equations. Each collocation point holds mechanical–thermal variables. The high order algorithm and the homotopy transformation allow reducing the number of tangent matrices to decompose and to avoid iterative procedure. Comparisons with the classical iterative solver (Jamal et al. in J Comput Mech 28:375–380, 2002) are performed. Numerical results reveal that a few number of matrix factorization is needed with the proposed approach, decreasing the computation time.