We analyze the far-infrared (FIR) properties of ~5000 star-forming galaxies at z<4.5, drawn from the deepest, super-deblended catalogs in the GOODS-N and COSMOS fields. We develop a novel panchromatic spectral energy distribution fitting algorithm, Stardust, that models the emission from stars, active galactic nuclei (AGNs), and infrared dust emission, without relying on energy balance assumptions. Our code provides robust estimates of the UV- optical and FIR physical parameters, such as the stellar mass (M*), dust mass (Mdust), infrared luminosities (LIR) arising from AGN and star formation activity, and the average intensity of the interstellar radiation field (<U>). Through a set of simulations we quantify the completeness of our data in terms of Mdust, LIR, and <U> and subsequently characterize the distribution and evolution of these parameters with redshift. We focus on the dust-to-stellar mass ratio (fdust), which we parameterize as a function of cosmic age, stellar mass, and specific star formation rate. The fdust is found to increase by a factor of 10 from z=0 to z=2 and appears to remain flat at higher z, mirroring the evolution of the gas fraction. We also find a growing fraction of warm to cold dust with increasing distance from the main sequence, indicative of more intense interstellar radiation fields, higher star formation efficiencies, and more compact star-forming regions for starburst galaxies. Finally, we construct the dust mass functions (DMFs) of star-forming galaxies up to z=1 by transforming the stellar mass function to DMF through the scaling relations derived here. The evolution of fdust and the recovered DMFs are in good agreement with the theoretical predictions of the Horizon-AGN and IllustrisTNG simulations.