  We present a finite volume method preserving the invariant region property
(IRP) for the reaction-diffusion systems with quasimonotone functions,
including nondecreasing, decreasing, and mixed quasimonotone systems. The
diffusion terms and time derivatives are discretized by a finite volume method
satisfying the discrete maximum principle (DMP) and the backward Euler method,
respectively. The discretization leads to an implicit and nonlinear scheme, and
it is proved to preserve the invariant region property unconditionally. We
construct an iterative algorithm and prove the invariant region property ar
each iteration step. Numerical examples are shown to confirm the accuracy and
invariant region property of our scheme.
