We consider the problem of heat transport by vibrational modes between Langevin thermostats connected by a central device. The latter is anharmonic and can be subject to large temperature difference and thus be out of equilibrium. We develop a classical formalism based on the equation of motion method, the fluctuation–dissipation theorem and the Novikov theorem to describe heat flow in a multi-terminal geometry. We show that it is imperative to include a quartic term in the potential energy to insure stability and to properly describe thermal expansion. The latter also contributes to leading order in the thermal resistance, while the usually adopted cubic term appears in the second order. This formalism paves the way for accurate modeling of thermal transport across interfaces in highly non-equilibrium situations beyond perturbation theory.