A procedure for fitting multi-exponential functions to experimental data is described. It is fast, requires no initial parameter estimates and is particularly suited to sums of several closely spaced exponentials. The method comprises the application of three well tried numerical techniques: (i) the signal is smoothed by representing it as an abbreviated Legendre series; (ii) the coefficients of a certain kind of differential equation are determined such that it's solution is the closest fit to the smoothed signal; and (iii) the amplitudes of the exponential components are determined, given the calculated values of the exponential rate constants. The method is computationally efficient, since determination of amplitudes and exponents involves the use of linear techniques, and therefore does not require multiple iterations, and the smoothed signal is contained in a handful of coefficients rather than as a lengthy time series. The severe ill-conditioning that is unavoidable in this problem is contained within the well-understood procedures of inverting a matrix and determining the roots of a polynomial. This method is particularly appropriate for analysis of data that may be modelled by a scheme of linked first-order reactions, describing for example the stochastic behaviour of ion channels, a chemical reaction, or the uptake and distribution of a drug within body compartments.