  We model the evolution of IR galaxies using a phenomenological approach to
match the observed source counts at different IR wavelengths. We introduce a
new algorithm for reproducing source counts based on direct integration of
probability distributions rather than Monte-Carlo sampling. We construct a
simple model for the evolution of the luminosity function and the colour
distribution of IR galaxies which utilizes a minimum number of free parameters.
Moreover we analyze how each of these parameters is constrained by
observational data. The model is based on pure luminosity evolution and adopts
the Dale & Helou SED templates. We find that the 850um source counts and their
redshift distribution depend strongly on the shape of the luminosity evolution
function, but only weakly on the details of the SEDs. We derive the best-fit
evolutionary model using the 850um counts and redshift distribution as
constraints. Moreover our best-fit shows a flattening of the faint end of the
luminosity function towards high redshifts and requires a colour evolution
which implies the typical dust temperatures of objects with the same
luminosities to decrease with redshift. We compare our best-fit model to
observed source counts at shorter and longer wavelengths which indicates our
model reproduces the 70um and 1100um source counts remarkably well, but
under-produces the counts at intermediate wavelengths. Analysis reveals that
the discrepancy arises at low redshifts, indicating that revision of the
adopted SED library towards lower dust temperatures (at a fixed infrared
luminosity) is required. This modification is equivalent to a population of
cold galaxies existing at low redshifts, as also indicated by recent Herschel
results, which are underrepresented in IRAS sample. We show that the modified
model successfully reproduces the source counts in a wide range of IR and submm
wavelengths.
