# 2D reaction-diffusion system simulation

I am a complete beginner in numerical simulation and I am pretty lost about how to tackle this problem.

I have been trying for some time to find the steady state (or simulate), the following system

$$\dot S = D_S \Delta S + \rho_S \frac{ F^n } { F^n + 1} - \delta_S S -\delta_C P^2 S$$

$$\dot P = D_P \Delta P + \rho_P - \delta_P P + (v - 2 \delta_C) P ^ 2 S$$

$$\dot F = \Delta + \rho_F \frac{ 1 } { (P^2S)^n + 1 } - \delta_F F$$

on a mesh that looks like this (with no boundary conditions given):

I realize that this systems contains the temporal derivatives, so I'm not explicitly going for the steady state, but when I took out the TransientTerm it complained that it was Exactly singular.

S = CellVariable(mesh = mesh, hasOld=True)
P = CellVariable(mesh = mesh, hasOld=True)
F = CellVariable(mesh = mesh, hasOld=True)

sdot = TransientTerm(var=S) == DiffusionTerm(coeff=Ds, var=S) + \
rs * pow(F,2) / (pow(F,2) + 1) - ds * S - dc * pow(P, 2) * S

pdot = TransientTerm(var=P) == DiffusionTerm(coeff=Dp, var=P) + \
rp  - dp * P  + (v - 2. * dc) * S * pow(P, 2)

fdot = TransientTerm(var=F) == DiffusionTerm(coeff=1, var=F) + \
rf / (pow(S * pow(P, 2), 2)  + 1.0 ) - df * F

dt = 0.005
coupled.sweep(dt = dt)

for i in range(0, 100000):

v1.plot('plots/SHH_'+str(dt * i)+'.png')
v2.plot('plots/Ptc_' + str(dt * i)+'.png')
v3.plot('plots/FGF_'+str(dt * i)+'.png')

S.updateOld()
P.updateOld()
F.updateOld()

dt = 0.005

for i in range(0, 100000):

v1.plot('plots/SHH_'+str(dt * i)+'.png')
v2.plot('plots/Ptc_' + str(dt * i)+'.png')
v3.plot('plots/FGF_'+str(dt * i)+'.png')

S.updateOld()
P.updateOld()
F.updateOld()
coupled.sweep(dt = dt)


The simulation initially appears to be proceeding as normal, but eventually, regardless of the size of my time-step, the $S$ values and $P$ values drop below zero and subsequently $S$ blows up do its large production rate $\rho_s = 1600$, causing an overflow error.

I should mention that when I simulated this problem in 1D using pdepe, the simulation was find and converged to a steady state.

How do I debug this, and is there other problems with my simulation? Should I try a different package such as Firedrake? If so, how do I solve the convergence problems there?