This article constructs and analyzes a numerical method for a time-dependent weakly coupled system of two singularly perturbed semilinear reaction–diffusion equations. The source terms in both equations have discontinuities in the spatial variables along the interface x=d,d∈Ω≔(0,1). The highest order spatial derivatives in the first and second equations are multiplied by positive perturbation parameters ɛ1 and ɛ2, respectively, which could be arbitrarily small. In the solution of the considered problem, boundary and interior layers appear in narrow regions of the domain. According to the boundary and interior layers, the solution is decomposed into regular and singular components, and precise bounds on the solution and its derivatives are given. The domain is discretized using an appropriate Shishkin mesh. For the mesh points which are not on the interface, the problem is discretized using a central difference method in space and backward Euler in time; for the mesh points which are on the interface, a special finite difference scheme is constructed. Parameters-uniform ((ɛ1,ɛ2)-uniform) error estimates in “maximum norm” have been obtained. It is proved that the method is parameters-uniformly convergent of first-order in time and almost second-order in space concerning perturbation parameters. Numerical experiments are conducted to demonstrate the efficiency of the method.