  In this article we present robust, efficient and accurate fully implicit
time-stepping schemes and nonlinear solvers for systems of reaction-diffusion
equations. The applications of reaction-diffusion systems is abundant in the
literature, from modelling pattern formation in developmental biology to cancer
research, wound healing, tissue and bone regeneration and cell motility.
Therefore, it is crucial that modellers, analysts and biologists are able to
solve accurately and efficiently systems of highly nonlinear parabolic partial
differential equations on complex stationary and sometimes continuously
evolving domains and surfaces. The main contribution of our paper is the study
of fully implicit schemes by use of the Newton method and the Picard iteration
applied to the backward Euler, the Crank-Nicolson (and its modifications) and
the fractional-step theta methods. Our results conclude that the
fractional-step theta method coupled with a single Newton iteration at each
timestep is as accurate as the fully adaptive Newton method; and both
outperform the Picard iteration. In particular, the results strongly support
the observation that a single Newton iteration is sufficient to yield as
accurate results as those obtained by use of an adaptive Newton method. This is
particularly advantageous when solving highly complex nonlinear partial
differential equations on evolving domains and surfaces. To validate our
theoretical results, various appropriate numerical experiments are exhibited on
stationary planary domains and in the bulk of stationary surfaces.
