  In this paper we present a fast algorithm for the numerical solution of
systems of reaction-diffusion equations, $\partial_t u + a \cdot \nabla u =
\Delta u + F (x, t, u)$, $x \in \Omega \subset \mathbf{R}^3$, $t > 0$. Here,
$u$ is a vector-valued function, $u \equiv u(x, t) \in \mathbf{R}^m$, $m$ is
large, and the corresponding system of ODEs, $\partial_t u = F(x, t, u)$, is
stiff. Typical examples arise in air pollution studies, where $a$ is the given
wind field and the nonlinear function $F$ models the atmospheric chemistry.
