The high mobility, density and multi-path evident in modern wireless systems makes the channel highly non-stationary. This causes temporal variation in the channel distribution that leads to the existence of time-varying joint interference across multiple degrees of freedom (DoF, e.g., users, antennas, frequency and symbols), which renders conventional precoding sub-optimal in practice. In this work, we derive a High-Order Generalization of Mercer's Theorem (HOGMT), which decomposes the multi-user non-stationary channel into two (dual) sets of jointly orthogonal subchannels (eigenfunctions), that result in the other set when one set is transmitted through the channel. This duality and joint orthogonality of eigenfuntions ensure transmission over independently flat-fading subchannels. Consequently, transmitting these eigenfunctions with optimally derived coefficients eventually mitigates any interference across its degrees of freedoms and forms the foundation of the proposed joint spatio-temporal precoding. The transferred dual eigenfuntions and coefficients directly reconstruct the data symbols at the receiver upon demodulation, thereby significantly reducing its computational burden, by alleviating the need for any complementary post-coding. Additionally, the eigenfunctions decomposed from the time-frequency delay-Doppler channel kernel are paramount to extracting the second-order channel statistics, and therefore completely characterize the underlying channel. We evaluate this using a realistic non-stationary channel framework built in Matlab and show that our precoding achieves ${\geqslant}$4 orders of reduction in BER at SNR${\geqslant}15$dB in OFDM systems for higher-order modulations and less complexity compared to the state-of-the-art precoding.