The aim of this paper is twofold: the first aim is to formulate and validate a multi-scale discrete Boltzmann method (DBM) based on density functional kinetic theory for thermal multiphase flow systems, ranging from continuum to transition flow regime; the second aim is to present some new insights into the thermo-hydrodynamic non-equilibrium (THNE) effects in the phase separation process. Methodologically, for bulk flow, DBM includes three main pillars: (i) the determination of the fewest kinetic moment relations, which are required by the description of significant THNE effects beyond the realm of continuum fluid mechanics; (ii) the construction of an appropriate discrete equilibrium distribution function recovering all the desired kinetic moments; (iii) the detection, description, presentation and analysis of THNE based on the moments of the non-equilibrium distribution ( $f-f^{(eq)}$). The incorporation of appropriate additional higher-order thermodynamic kinetic moments considerably extends the DBM's capability of handling larger values of the liquid–vapour density ratio, curbing spurious currents, and ensuring mass/momentum/energy conservation. Compared with the DBM with only first-order THNE (Gan et al., Soft Matt., vol. 11 (26), 2015, pp. 5336–5345), the model retrieves kinetic moments beyond the third-order super-Burnett level, and is accurate for weak, moderate and strong THNE cases even when the local Knudsen number exceeds $1/3$. Physically, the ending point of the linear relation between THNE and the concerned physical parameter provides a distinct criterion to identify whether the system is near or far from equilibrium. Besides, the surface tension suppresses the local THNE around the interface, but expands the THNE range and strengthens the THNE intensity away from the interface through interface smoothing and widening.