We have adapted a simple method of numerical integration to predict the equilibrium state of a population of components undergoing reversible association according to the Law of Mass Action. Its particular application is to populations of protein molecules in aqueous solution. The method is based on Euler integration but employs an adaptive step size: the time increment being reduced if it would make the concentration of any component negative and increased while the concentration of any component changes at greater than a specified rate. Parameters of the algorithm have been optimized empirically using a model set of binding equilibria with dissociation constants ranging from 10(-5) M to 10(-9) M. The method obtains the solution to a set of binding equilibria more rapidly than the conventional initial value methods (simple Euler, 4th order Runge-Kutta and variable-step Runge-Kutta methods were tested) for the same accuracy. A computer code in standard C is presented.