  A new class of implicit-explicit (IMEX) methods combined with a p-adaptive
mixed finite element formulation is proposed to simulate the diffusion of
reacting species. Hierarchical polynomial functions are used to construct an
$H(\mathrm{Div})$-conforming base for the flux vectors, and a non-conforming
$L^2$ base for the mass concentration of the species. The mixed formulation
captures the distinct nonlinearities associated with the constitutive flux
equations and the reaction terms. The IMEX method conveniently treats these two
sources of nonlinearity implicitly and explicitly, respectively, within a
single time-stepping framework. The combination of the p-adaptive mixed
formulation and the IMEX method delivers a robust and efficient algorithm. The
proposed methods eliminate the coupled effect of mesh size and time step on the
algorithmic stability. A residual based a posteriori error estimate that
provides an upper bound of the natural error norm is derived. The availability
of such estimate which can be obtained with minimal computational effort and
the hierarchical construction of the finite element spaces allow for the
formulation of an efficient p-adaptive algorithm. A series of numerical
examples demonstrate the performance of the approach. It is shown that the
method with the p-adaptive strategy accurately solves problems involving
travelling waves, and those with discontinuities and singularities. The
flexibility of the formulation is also illustrated via selected applications in
pattern formation and electrophysiology.
