  We develop numerical methods for reaction-diffusion systems based on the
equations of fluctuating hydrodynamics (FHD). While the FHD formulation is
formally described by stochastic partial differential equations (SPDEs), it
becomes similar to the reaction-diffusion master equation (RDME) description
when those SPDEs are spatially discretized and reactions are modeled as a
source term having Poisson fluctuations. However, unlike the RDME, the FHD
description naturally extends from the regime where fluctuations are strong,
i.e., each hydrodynamic cell has few (reactive) molecules, to regimes with
moderate or weak fluctuations, and ultimately to the deterministic limit. By
treating diffusion implicitly, we avoid the severe restriction on time step
size that limits all methods based on explicit treatments of diffusion, and
construct numerical methods that are more efficient than RDME methods, without
compromising accuracy. We find that an implicit midpoint tau leaping scheme
attains second-order weak accuracy for the linearized FHD, and gives an
accurate and stable structure factor for a time step size an order of magnitude
larger than the hopping time scale of diffusing molecules. We study the
numerical accuracy of our methods for the Schl\"ogl reaction-diffusion model
both in and out of thermodynamic equilibrium. We demonstrate and quantify the
importance of thermodynamic fluctuations to the formation of a two-dimensional
Turing-like pattern, and examine the effect of fluctuations on
three-dimensional chemical front propagation. By comparing FHD simulations to
deterministic reaction-diffusion simulations, we show that fluctuations
accelerate pattern formation in spatially homogeneous systems, and lead to a
qualitatively-different disordered pattern behind a traveling wave.
