The parabolic differential equations with Neumann conditions and history functions have the following general form :
$$\left\{\begin{array}{lc} \dfrac{\partial u(t,x)}{\partial t}=d_1\triangle u(t,x)+u(t,x)\left(r_1-a_{11}u(t-\tau,x)-a_{12}v(t,x)\right),& t>0,x\in\Omega \\ \dfrac{\partial v(t,x)}{\partial t}=d_2\triangle v(t,x)+v(t,x)\left(-r_2+a_{21}u(t,x)-a_{22}v(t,x)\right),& t>0,x\in\Omega\\ \dfrac{\partial u}{\partial n}=\dfrac{\partial v}{\partial n}=0,\qquad t\ge0,x\in\partial\Omega \quad\text{(Neumann boundary conditions)}\\ u(t,x)=\phi(t,x)\ge 0,\qquad v(t,x)=\psi(t,x)\ge 0, &(t,x)\in[-\tau,0]\times\Omega \end{array} \right.$$
Where $d_1,d_2,r_1,r_2,a_{ij}$ are constants, $\Omega$ can be interval like $[0,100]$, $\psi$ and $\phi$ are known functions, e.g., $x\cdot \sin^4(x)$, $9\cos^22x$
Matlab, Python or R, only supports ordinary differential equations with time delay. How to solve such delayed partial differential equations using the already available tools in these mathematical softwares?
To my knowledge, in Python, the best tool, considering documentation, code quality, and the fact that it is being actively maintained is JitCDDE.
It's solver uses just-in-time compilation and Sympy expressions to define the model, so that algebraic optimizations to the formula are done automatically before the code is compiled leading to very efficient computations.