I debug the source code
R = 3.95; eps = 10^-3;
d2 = 0.03;
d3 = 11;
alpha1 = 0.2;
alpha2 = 0.12/60;
alpha3 = 1;
beta1 = 0.266;
beta2 = 0.28;
beta3 = 1;
gamma1 = 0.2667;
gamma2 = 0.35;
delta1 = 0.00297;
delta2 = 0.35;
eq1 = {\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(mT[r, theta, phi,
t]\)\) == (alpha1*mBG[r, theta, phi, t] + alpha2)*
mD[r, theta, phi, t] - alpha3*mT[r, theta, phi, t] +
beta1*mBG[r, theta, phi, t]*cD[r, theta, phi, t] +
d2*(1./(R^2*(Sin[phi])^2) \!\(
\*SubscriptBox[\(\[PartialD]\), \(theta, theta\)]\(mT[r, theta, phi,
t]\)\) + 1./(R^2*Sin[phi])*Cos[phi]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi\)]\(mT[r, theta, phi, t]\)\) +
1./R^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi, phi\)]\(mT[r, theta, phi,
t]\)\)),
\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(mD[r, theta, phi,
t]\)\) == -(alpha1*mBG[r, theta, phi, t] + alpha2)*
mD[r, theta, phi, t] + alpha3*mT[r, theta, phi, t] +
beta2*cD[r, theta, phi, t] - beta3*mD[r, theta, phi, t] +
d2*(1./(R^2*(Sin[phi])^2) \!\(
\*SubscriptBox[\(\[PartialD]\), \(theta, theta\)]\(mD[r, theta, phi,
t]\)\) + 1./(R^2*Sin[phi])*Cos[phi]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi\)]\(mD[r, theta, phi, t]\)\) +
1./R^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi, phi\)]\(mD[r, theta, phi,
t]\)\)),
\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(mB[r, theta, phi, t]\)\) ==
gamma1*mT[r, theta, phi, t]*cB[r, theta, phi, t] -
gamma2*mB[r, theta, phi, t] -
delta1*mB[r, theta, phi, t]*cG[r, theta, phi, t] +
delta2*mBG[r, theta, phi, t] +
d2*(1./(R^2*(Sin[phi])^2) \!\(
\*SubscriptBox[\(\[PartialD]\), \(theta, theta\)]\(mB[r, theta, phi,
t]\)\) + 1./(R^2*Sin[phi])*Cos[phi]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi\)]\(mB[r, theta, phi, t]\)\) +
1./R^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi, phi\)]\(mB[r, theta, phi,
t]\)\)),
\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(mBG[r, theta, phi, t]\)\) ==
delta1*mB[r, theta, phi, t]*cG[r, theta, phi, t] -
delta2*mBG[r, theta, phi, t] +
d2*(1./(R^2*(Sin[phi])^2) \!\(
\*SubscriptBox[\(\[PartialD]\), \(theta, theta\)]\(mBG[r, theta, phi,
t]\)\) + 1./(R^2*Sin[phi])*Cos[phi]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi\)]\(mBG[r, theta, phi, t]\)\) +
1./R^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi, phi\)]\(mBG[r, theta, phi,
t]\)\))};
eq2 = {\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(cD[r, theta, phi, t]\)\) ==
d3*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(r, r\)]\(cD[r, theta, phi, t]\)\) +
2./r*\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(cD[r, theta, phi, t]\)\) +
1./(r^2*(Sin[phi])^2) \!\(
\*SubscriptBox[\(\[PartialD]\), \(theta, theta\)]\(cD[r, theta, phi,
t]\)\) + 1./(r^2*Sin[phi])*Cos[phi]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi\)]\(cD[r, theta, phi, t]\)\) +
1./r^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi, phi\)]\(cD[r, theta, phi,
t]\)\)),
\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(cB[r, theta, phi, t]\)\) ==
d3*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(r, r\)]\(cB[r, theta, phi, t]\)\) +
2./r*\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(cB[r, theta, phi, t]\)\) +
1./(r^2*(Sin[phi])^2) \!\(
\*SubscriptBox[\(\[PartialD]\), \(theta, theta\)]\(cB[r, theta, phi,
t]\)\) + 1./(r^2*Sin[phi])*Cos[phi]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi\)]\(cB[r, theta, phi, t]\)\) +
1./r^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi, phi\)]\(cB[r, theta, phi,
t]\)\)),
\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(cG[r, theta, phi, t]\)\) ==
d3*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(r, r\)]\(cG[r, theta, phi, t]\)\) +
2./r*\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(cG[r, theta, phi, t]\)\) +
1./(r^2*(Sin[phi])^2) \!\(
\*SubscriptBox[\(\[PartialD]\), \(theta, theta\)]\(cG[r, theta, phi,
t]\)\) + 1./(r^2*Sin[phi])*Cos[phi]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi\)]\(cG[r, theta, phi, t]\)\) +
1./r^2*\!\(
\*SubscriptBox[\(\[PartialD]\), \(phi, phi\)]\(cG[r, theta, phi,
t]\)\))
};
intis = {mT[r, theta, phi, 0] == 0, mD[r, theta, phi, 0] == 0,
mB[r, theta, phi, 0] == 0, mBG[r, theta, phi, 0] == 0,
cD[r, theta, phi, 0] == 0.001, cB[r, theta, phi, 0] == 0.001,
cG[r, theta, phi, 0] == 0.001};
bc = {mT[r, theta, eps, t] == 0.001,
mT[r, theta, 2*Pi, t] == 0.001,
mT[r, eps, phi, t] == 0.001,
mT[r, Pi, phi, t] == 0.001,
mD[r, theta, eps, t] == 0.001,
mD[r, theta, 2*Pi, t] == 0.001,
mD[r, eps, phi, t] == 0.001,
mD[r, Pi, phi, t] == 0.001,
mB[r, theta, eps, t] == 0.001,
mB[r, theta, 2*Pi, t] == 0.001,
mB[r, eps, phi, t] == 0.001,
mB[r, Pi, phi, t] == 0.001,
mBG[r, theta, eps, t] == 0.001,
mBG[r, theta, 2*Pi, t] == 0.001,
mBG[r, eps, phi, t] == 0.001,
mBG[r, Pi, phi, t] == 0.001,
cD[r, theta, eps, t] == 0.001,
cD[r, theta, 2*Pi, t] == 0.001,
cD[r, eps, phi, t] == 0.001,
cD[r, Pi, phi, t] == 0.001,
cD[R, theta, phi, t] == 0.001,
d3*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(cD[r, theta, phi, t]\)\) /.
r -> R) == (-(beta1*mBG[R, theta, phi, t] + beta2)*
cD[R, theta, phi, t] + beta3*mD[R, theta, phi, t]),
cB[r, theta, eps, t] == 0.001,
cB[r, theta, 2*Pi, t] == 0.001,
cB[r, eps, phi, t] == 0.001,
cB[r, Pi, phi, t] == 0.001,
cB[R, theta, phi, t] == 0.001,
d3*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(cB[r, theta, phi, t]\)\) /.
r -> R) == (-gamma1*mT[R, theta, phi, t]*
cB[R, theta, phi, t] +
gamma2*mB[R, theta, phi, t]),
cG[r, theta, eps, t] == 0.001,
cG[r, theta, 2*Pi, t] == 0.001,
cG[r, eps, phi, t] == 0.001,
cG[r, Pi, phi, t] == 0.001,
cG[R, theta, phi, t] == 0.001,
d3*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(cG[r, theta, phi, t]\)\) /.
r -> R) == (-delta1*mB[R, theta, phi, t]*
cG[R, theta, phi, t] +
delta2*mBG[R, theta, phi, t])
};
usol = NDSolveValue[
Flatten[{eq1, eq2, intis, bc}], {mT, mD, mB, mBG, cD, cB, cG}, {r,
eps, R}, {theta, eps, Pi}, {phi, eps, 2*Pi}, {t, 0, 0.01}]
var = {mT, mD, mB, mBG, cD, cB, cG};
With[{r = R, t = .01},
Table[Plot3D[
usol[[i]][R, theta, phi, .01], {theta, eps, Pi}, {phi, eps, 2*Pi},
PlotLabel -> var[[i]], PlotRange -> All, Mesh -> None,
ColorFunction -> Hue], {i, 1, Length[var]}]] // Quiet

y[x] /. x -> R == 0and observe the output, and think about why you obtain this. – xzczd Feb 18 '19 at 14:28(y[x] /.x->R) == 0. However, now I get the error:Transpose::nmtx: The first two levels of {mT,mD,mB,mBG,NDSolvexs$1044,NDSolvexs$1043,NDSolvexs$1045} cannot be transposed.Here it seems something is wrong with the functions aftermGB`, but I cannot figure it out. Why the name of the functions in this list involve NDSolve? – MOON Feb 18 '19 at 14:55Withlike in e.g. this post to simplify the code. – xzczd Feb 18 '19 at 15:21NDSolveinstead ofNDSolveValue. Maybe you can fix that first. -- Minor style remark: I think{f1, f2} = ...is a bad idea, because then what were variables ineqnsetc. now have values, making those previously defined things invalid. It makes the code less robust, hard to play with and debug, imo. – Michael E2 Feb 18 '19 at 16:41NDSolveandNDSolveValue? – MOON Feb 19 '19 at 14:48NDSolveis the critical one and perhaps causes the others messages inNDSolveValue. I'm suggesting you look at theNDSolvemessage, which indicates the problem is ill-posed. I don't know how to fix it, so I'm suggesting you might come up with an appropriate change. Alex Trounev basically guesses at a possible way to change the problem, which does in fact create a new problem thatNDSolvewill solve. But I don't know if the new problem is the one you want to solve or not. – Michael E2 Feb 19 '19 at 19:51