The eikonal approximation (instanton technique) is applied to the problem of large fluctuations of the number of species in spatially homogeneous chemical reactions with the probability density distribution described by a master equation. For both autocatalytic and nonautocatalytic reactions, the analysis of the distribution about a stable stationary state and of the transitions between coexisting stable states comes, to logarithmic accuracy, to the analysis of Hamiltonian dynamics of an auxiliary dynamical system. The latter can be done explicitly in a few cases, including one-species systems, systems with detailed balance, and systems close to the bifurcation points where the number of the stable states changes. In the last case, the fluctuations display universal features, and, for saddle-node bifurcation points, the logarithm of the probability of escape from the metastable state (per unit time) is proportional to the distance to the bifurcation point (in the parameter space) raised to the power 3/2. We compare the eikonal approximation for the stationary distribution of a master equation to Monte Carlo numerical solutions for two chemical two-variable systems with multiple stationary states, where none of the cited restrictions exists. For one of the systems in the pattern of optimal paths we observe caustics emanating from the saddle point.