High order parameter-robust numerical method for a system of (M>=2) coupled singularly perturbed parabolic reaction-diffusion problems

We present a high order parameter-robust numerical method for a system of (M>=2) coupled singularly perturbed parabolic reaction-diffusion problems. A small perturbation parameter {\epsilon} is multiplied with the second order spatial derivatives in all the equations. The parabolic boundary layer appears in the solution of the problem when the perturbation parameter {\epsilon} tends to zero. To obtain a high order approximation to the solution of this problem, we propose a numerical method that employs the Crank-Nicolson method on an uniform mesh in time direction, together with a hybrid finite difference scheme on a generalized Shishkin mesh in spatial direction. We prove that the resulting method is parameter-robust or \epsilon-uniform of second order in time and almost fourth order in spatial variable, if the discretization parameters satisfy a non-restrictive relation. Numerical experiments are presented to validate the theoretical results and also indicate that the relation between the discretization parameters is not necessary in practice.