Covariance functions are the equivalent of covariance matrices for traits with many, potentially infinitely many, records in which the covariances are defined as a function of age or time. They can be fitted for any source of variation, e.g. genetic, permanent environment or phenotypic. A suitable family of functions for covariance functions are orthogonal polynomials. These give the covariance between measurements at any two ages as a higher order polynomial of the ages at recording. Polynomials can be fitted to full or reduced order. The former is equivalent to a multivariate analysis estimating covariance components. A reduced order fit involves less parameters and smoothes out differences in estimates of covariances. It gives predicted covariance matrices of rank equal to the order of fit. The coefficients of covariance functions can be estimated by restricted maximum likelihood through a reparameterisation of existing algorithms to estimate covariance components. For a simple animal model with equal design matrices for all traits, computational requirements to estimate covariance functions are proportional to the order of fit for the genetic covariance function. Applications to simulated data and a set of beef cattle data are shown.