The method of least squares (LS) yields exact solutions for the adjustable parameters when the number of data values n equals the number of parameters "p". This holds also when the fit model consists of "m" different equations and "m = p", which means that LS algorithms can be used to obtain solutions to systems of equations. In particular, nonlinear LS solves systems of nonlinear equations. An important example in chemistry is the case of reagents whose concentrations are coupled through multiple equilibrium relations. The capability of nonlinear LS in this application is examined for three programming environments, Excel Solver, FORTRAN, and KaleidaGraph, on a number of equilibrium problems having up to 10 unknown concentrations. FORTRAN and KaleidaGraph perform well in all the examples, but Solver presents difficulties that render it inadequate in several cases unless the problem is reformulated in terms of a smaller number of adjustable concentrations. When the input quantities (equilibrium constants, prepared concentrations) have uncertainty, the calculations can also be used to propagate these uncertainties into the derived quantities.