NDSolve::bcedge: Boundary condition c[t,5]==Cout is not specified on a single edge of the boundary of the computational domain. >>

I'd like to plot $\frac{\partial}{\partial t}c=\frac{d}{r^2}\frac{\partial}{\partial r}(r^2 \frac{\partial}{\partial r}c) \equiv\Delta c $ with the initial condition $c(0,r)=c_{0}$ and the boundary conditions $\frac{\partial}{\partial r}c(t,0)=0$ and $c(t,R\in\mathbb{R})=c_{out}$ where $R$ is the radius of a circle. I know the analytical solution and I know how the profile looks. I'd like to use/learn Mathematica because it often helps if you can make a quick plot of unknown shapes.

NDSolve[{D[c[t, r], t] == d/(r^2) D[((r^2) D[c[t, r], r]), r],
Derivative[0, 1][c][t, 0] == 0, c[t, 5] == cout,c[0, r] == c0}, 
c, {t, 0, 10}, {r, -5, 5}]
  • $\begingroup$ Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$
    – bbgodfrey
    Apr 8, 2015 at 21:50

1 Answer 1


The specific error occurs because the range in r is given as {r, -5, 5} but the boundary condition Derivative[0, 1][c][t, 0] == 0 is given at r = 0, which is not a boundary. I imagine that {r, 0, 5} is meant, which eliminates the error. However, the equation is singular at r = 0, which creates other errors. This is a common issue in spherical coordinates. An easy work-around is to displace the inner radial boundary slightly, say to r = 0.01. Finally, all constants need to be specified. In all, I modified the code to

cout = 1; c0 = 2; d = 1;
ans = NDSolveValue[{D[c[t, r], t] == d/(r^2) D[((r^2) D[c[t, r], r]), r], 
   Derivative[0, 1][c][t, 0.01] == 0, c[t, 5] == cout, c[0, r] == c0},
   c, {t, 0, 10}, {r, 0.1, 5}];

and then plotted the solution

Plot3D[ans[t, r], {t, 0, 10}, {r, 0.01, 5}, AxesLabel -> {t, r, c}]

Mathematica graphics

This may get you started.

  • $\begingroup$ Thanks for the example code and explanation! Actually my version does not support NDSolveValue and NDSolve returns NDSolve::ibcinc: Warning: Boundary and initial conditions are inconsistent and the plot is empty. Your guess on the BC was right. My attempt by choosing {r,-5,5} was to get a symmetrically continued function plot over r<0 (mirroring at the c-t-plane). Thanks a lot for your welcoming comment. $\endgroup$
    – dkeck
    Apr 9, 2015 at 6:58
  • $\begingroup$ NDSolve works just as well, and you usually can ignore the warning. See the documentation of NDSolve on how to extract and plot the results it produces. $\endgroup$
    – bbgodfrey
    Apr 9, 2015 at 12:16
  • $\begingroup$ Specifically, replace NDSolveValue[... by First@NDSolve[..., and ans[t, r] by (c /. ans)[t, r] in Plot3D. $\endgroup$
    – bbgodfrey
    Apr 9, 2015 at 12:44
  • $\begingroup$ I already tried ans = NDSolve[{D[c[t, r], t] == d/(r^2) D[((r^2) D[c[t, r], r]), r], Derivative[0, 1][c][t, 0.01] == 0.01, c[t, 5] == cout, c[0, r] == c0}, c, {t, 0, 10}, {r, 0.01, 5}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 1000}}] with Plot3D[Evaluate[c[t, r] /. ans], {t, 0, 10}, {r, 0.01, 5}, AxesLabel -> {t, r, c}, PlotRange -> {1, 3}] In both cases I only see a 2D plane constant to 2. Where MethodOfLines option was used to overcome the inconsistent boundary warning. $\endgroup$
    – dkeck
    Apr 9, 2015 at 20:26
  • $\begingroup$ What values are you using for your constants?. If cout=c0=2, the answer is 2 everywhere. $\endgroup$
    – bbgodfrey
    Apr 10, 2015 at 1:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.