# Getting error when trying to solve two coupled differential equations with NDSolve

Mathematica newbie here. I have been trying to use NDSolve to solve two coupled differential equations. Each equation on its own (minus the coupling terms) solves fine, but when I try them both together I get an error. I apologize in advance if my question is stupid or this is not the place to ask it. Here is my code:

FIRST EQUATION:

xb = 1;
kcat = 1;
v = 1;
kp = 1;
kd = 1;
Dp = 1;
Dc = 1;
kform = 1;

NDSolve[
{D[P[x, t], t] == Dp*D[P[x, t], x, x] - kcat*(*C[x,t]*)P[x, t] +
(kp/2)*(1/Sqrt[π*v])*(Exp[-((x - xb)^2)/v] + Exp[-((x + xb)^2)/v]),
(D[P[x, t], x] /. x -> -xb) == 0,
(D[P[x, t], x] /. x -> xb) == 0,
P[x, 0] == 0},
P[x, t], {x, -xb, xb}, {t, 0, 50}]

Quit[]


SECOND EQUATION:

xb = 1;
kcat = 1;
v = 1;
kp = 1;
kd = 1;
Dp = 1;
Dc = 1;
kform = 1;

NDSolve[
{D[C[x, t], t] == Dc*D[C[x, t], x, x] - kd*C[x, t] +
kform*(1/Sqrt[π*v])*Exp[-(x^2)/v],
(D[C[x, t], x] /. x -> -xb) == 0,
(D[C[x, t], x] /. x -> xb) == 0,
C[x, 0] == 5},
C[x, t], {x, -xb, xb}, {t, 0, 50}]

Quit[]


BOTH:

xb = 1;
kcat = 1;
v = 1;
kp = 1;
kd = 1;
Dp = 1;
Dc = 1;
kform = 1;

NDSolve[
{D[P[x, t], t] == Dp*D[P[x, t], x, x] - kcat*C[x, t]*P[x, t],
D[C[x, t], t] == Dc*D[C[x, t], x, x],
(D[P[x, t], x] /. x -> -xb) == 0,
(D[P[x, t], x] /. x -> xb) == 0,
P[x, 0] == 0,
(D[C[x, t], x] /. x -> -xb) == 0,
(D[C[x, t], x] /. x -> xb) == 0,
C[x, 0] == 5},
{P[x, t], C[x, t]}, {x, -xb, xb}, {t, 0, 50}]

Quit[]


ERROR MESSAGE:

NDSolve::ndode: The equations are not differential equations or initial conditions in the dependent variables {P}.

• Since C is a built-in function, I would not recommend using it as the name of the function involved in your PDE. – J. M. is in limbo Feb 18 '16 at 2:59

I rewrote your equations with Q in place of C which eliminated the error message and made them reality solvable, but I don't think the solutions are much good. It is my belief that your boundary conditions are ill- formulated. Since I have no idea what system you are modeling, there is no way for me help you find correct boundary conditions. Indeed, for lack of domain knowledge, there is a good chance I could not help you even you were to describe the physical system behind your equations.

I will say that I look upon the boundary conditions given by

(D[P[x, t], x] /. x -> -xb) == 0
(D[P[x, t], x] /. x -> xb) == 0


with a great deal of suspicion.

So you can see what I found and perhaps draw some conclusions from it, I present my work.

xb = 1;
kcat = 1;
v = 1;
kp = 1;
kd = 1;
Dp = 1;
Dc = 1;
kform = 1;

p =
NDSolve[
{D[P[x, t], t] == Dp*D[P[x, t], x, x] - kcat*P[x, t] +
(kp/2)*(1/Sqrt[π*v])*(Exp[-((x - xb)^2)/v] + Exp[-((x + xb)^2)/v]),
(D[P[x, t], x] /. x -> -xb) == 0,
(D[P[x, t], x] /. x -> xb) == 0,
P[x, 0] == 0},
P, {x, -xb, xb}, {t, 0, 50}][[1, 1, 2]];

Plot3D[p[x, t], {x, -1., 1.}, {t, 0., 50.}, PlotRange -> All] q =
NDSolve[
{D[Q[x, t], t] == Dc*D[Q[x, t], x, x] - kd*Q[x, t] +
kform*(1/Sqrt[π*v])*Exp[-(x^2)/v],
(D[Q[x, t], x] /. x -> -xb) == 0,
(D[Q[x, t], x] /. x -> xb) == 0,
Q[x, 0] == 5},
Q, {x, -xb, xb}, {t, 0, 50}][[1, 1, 2]];

Plot3D[q[x, t], {x, -1., 1.}, {t, 0., 50.}, PlotRange -> All] {p2, q2} =
NDSolve[
{D[P[x, t], t] == Dp*D[P[x, t], x, x] - kcat*Q[x, t]*P[x, t],
D[Q[x, t], t] == Dc*D[Q[x, t], x, x],
(D[P[x, t], x] /. x -> -xb) == 0,
(D[P[x, t], x] /. x -> xb) == 0,
P[x, 0] == 0,
(D[Q[x, t], x] /. x -> -xb) == 0,
(D[Q[x, t], x] /. x -> xb) == 0, Q[x, 0] == 5},
{P, Q}, {x, -xb, xb}, {t, 0, 50}][[1, All, 2]];

Plot3D[p2[x, t], {x, -1., 1.}, {t, 0., 50.}, PlotRange -> All] Plot3D[q2[x, t], {x, -1., 1.}, {t, 0., 50.}, PlotRange -> All] 