In this work we present the bulk-surface finite element method (BSFEM) for solving coupled systems of
bulk-surface reaction–diffusion equations (BSRDEs) on stationary volumes. Such systems of coupled
bulk-surface partial differential equations arise naturally in biological applications and fluid dynamics,
for example, in modelling of cellular dynamics in cell motility and transport and diffusion of surfactants
in two phase flows. In this proposed framework, we define the surface triangulation as a collection of the
faces of the elements of the bulk triangulation whose vertices lie on the surface. This implies that the
surface triangulation is the trace of the bulk triangulation. As a result, we construct two finite element
spaces for the interior and surface respectively. To discretise in space we use piecewise bilinear elements
and the implicit second order fractional-step θ scheme is employed to discretise in time. Furthermore,
we use the Newton method to treat the nonlinearities. The BSFEM applied to a coupled system of BSRDEs
reveals interesting patterning behaviour. For a set of appropriate model parameter values, the surface
reaction–diffusion system is not able to generate patterns everywhere in the bulk except for a small
region close to the surface while the bulk reaction–diffusion system is able to induce patterning almost
everywhere. Numerical experiments are presented to reveal such patterning processes associated with
reaction–diffusion theory.