# Use FEM Solver for System of Nonlinear PDEs with Nonlinear Coefficients - Reaction-Diffusion

I am trying to solve a complex transient reaction-diffusion problem numerically in 2D. Let's take a simple system to start, say

$$A+B\rightarrow C \\ C \rightarrow A+B$$

with reaction rate $k_1$ and $k_2$, respectively.

Looking at a radially symmetric system $(r,z)$ with no flux boundary conditions, $[A]$ being generated within the domain with arbitrary spatial distribution, and fixed initial concentrations, I put together the following in Mathematica 10.1:

(* Geometric Params *)
r1 = 0; r2 = 50; z1 = 0; z2 = 0.2; rb = 1;

(* Diffusion Coefficients *)
dA = 1; dB = 1; dC = 1;

(* Rate Constants *)
k1 = 1; k2 = 2;

(* Generation Rate *)
generation = 10;

(* Helper Function for Operator *)
operator[f_, d_] := D[f, t] - d Laplacian[f, {r, z}];

sol = NDSolveValue[
{
operator[cA[t, r, z], dA] == -k1 cA[t, r, z] cB[t, r, z] +
k2 cC[t, r, z] + generation (1 - 1/(1 + Exp[-20 r + 20]))
,operator[cB[t, r, z], dB] == -k1 cA[t, r, z] cB[t, r, z] +
k2 cC[t, r, z]
,operator[cC[t, r, z], dC] ==
k1 cA[t, r, z] cB[t, r, z] - k2 cC[t, r, z]

,DirichletCondition[cA[t, r, z] == 0, True]
,DirichletCondition[cB[t, r, z] == 0, True]
,DirichletCondition[cC[t, r, z] == 0, True]

,cA[0, r, z] == 0
,cB[0, r, z] == 100
,cC[0, r, z] == 0

}, {cA, cB, cC}, {t, 0, 100}, {r, r1, r2}, {z, z1, z2}
, Method -> {"MethodOfLines", "TemporalVariable" -> t}
]


This gives me an error stating that "Nonlinear coefficients are not supported in this version of NDSolve".

Is there a work around to solve a system like this?

I tried playing around with some of the low level FEM solver functions, but get similar errors along the way.

I really appreciate the help!

Replace

,DirichletCondition[cA[t, r, z] == 0, True]
,DirichletCondition[cB[t, r, z] == 0, True]
,DirichletCondition[cC[t, r, z] == 0, True]


by

, DirichletCondition[cA[t, r, z] == 0, z == z1 || z == z2]
, DirichletCondition[cB[t, r, z] == 0, z == z1 || z == z2]
, DirichletCondition[cC[t, r, z] == 0, z == z1 || z == z2]
, DirichletCondition[cA[t, r, z] == 0, r == r1 || r == r2]
, DirichletCondition[cB[t, r, z] == 0, r == r1 || r == r2]
, DirichletCondition[cC[t, r, z] == 0, r == r1 || r == r2]


in order to specify boundary conditions in r and z, but not in t, which already are specified by cA[0, r, z] == 0, etc.

Even these boundary conditions generate a warning message,

NDSolveValue::ibcinc: Warning: boundary and initial conditions are inconsistent. >>


because cB[0, r, z] == 100 intersects other boundaries where cB ==0, but I do not believe this to be a serious issue. A sample plot for cB is

Plot3D[sol[[2]][t, r, .1], {t, 0, 100}, {r, r1, r2}, PlotRange -> All]


sol[[1]] and sol[[3]] are zero.