There are considered: a bounded domain \(\Omega\subset \mathbb{R}^p\) \((p= 1,2,\dots)\), a finite interval \([0, T]\), the sets \(D_T= (0, T]\times \Omega\), \(S_t= (0, T]\times \partial \Omega\), \(\overline D_T= [0, T]\times \overline\Omega\), \(J_i= [-r_i, 0]\) (\(r_i\), \(i= 1,\dots, n\), represent the finite time delays), \(Q^{(i)}_T= [- r_i, T]\times \overline\Omega\), \(Q_T= Q^{(1)}_T\times\cdots\times Q^{(n)}_T\), and the system \[ \partial u_i/\partial t- L_i u_i= f_i(t, x, u(t, x), u_t(t, x)),\quad i= 1,\dots, n,\quad \text{in } D_t,\tag{1} \] where \(u(t, x)\equiv (u_1(t, x),\dots, u_n(t, x))\), \(u_t(t, x)= (u_1(t- r_1, x),\dots, u_n(t- r_n, x))\) (\(u_i(t, x)\) is a density function), and for each \(i\), \(L_i\) is a uniformly elliptic operator given in the form \[ L_i u_i\equiv \sum^p_{j, k= 1} a^{(i)}_{j,k}(t, x) \partial^2 u_i/\partial x_j \partial x_k + \sum^p_{j= 1} b^{(i)}_j (t, x) \partial u_i/\partial x_j. \] Some boundary and initial conditions are also given on \(S_T\) and \(J_i\times \Omega\), \(i= 1,\dots, n\). The functions \(f_i(t, x, u, v)\) are Hölder continuous in \((t, x)\) and locally Lipschitz continuous in \((u, v)\). In addition to the system (1) the system \[ \partial u_i/\partial t- L_i u_i= f_i(t, x, u(t, x), u_t(t, x))+ \int_\Omega g_i(t, x, x', u(t, x'), u_t(t, x')) dx',\quad i= 1,\dots, n,\tag{2} \] is also considered, where \(g_i\) satisfy the same properties as \(f_i\). A monotone iterative scheme, using upper-lower solutions, is studied and some existence-comparison theorems for (1) and (2) are given. Some applications to several model problems arising from ecology and nuclear engineering are presented.