We present a Cartesian grid finite difference numerical method for solving a system of reaction-diffusion initial boundary value problems with Neumann type boundary conditions. The method utilizes adaptive time-stepping, which guarantees stability and non-negativity of the solutions. The latter property is critical for models in biology where solutions rep- resent physical measurements such as concentration. The level set representation of the boundary enables us to handle domains with complicated geometry with ease. We pro- vide numerical validation of our method on synthetic and biological examples. Empirical tests demonstrate second order convergence rate in the L1- and L2-norms, as well as in the L(infinity)-norm for many cases.