In this talk, I will discuss the two-loop calculation of the bare next-to-leading power soft functions appearing in the Drell-Yan process in the threshold kinematics. These soft functions are the last missing ingredient needed to verify the next-to-leading power Drell-Yan factorization formula to NNLO. I will briefly review the derivation of the bare factorization theorem and introduce the relevant basis of the power suppressed soft functions. The talk will be dedicated to the discussion of the technical aspects of the calculation. In particular, I will focus on the evaluation of the master integrals using the differential equations method.