# plot multiple curves : solution obtained from NDSolve

I want to plot the solution obtained from NDSolve as a multicurve plot like shown in the following sample image (image ref). Here, x1, x2 , x3, x4 .. are the solution computed at includePoints = {{10}, {20}, {30}, {40}, {50}} in my case. In the following code (ref):

Needs["NDSolveFEM"]
region = Line[{{0}, {100}}];
includePoints = {{10}, {20}, {30}, {40}, {50}};
mesh = ToElementMesh[region, "IncludePoints" -> includePoints,
"MaxCellMeasure" -> 0.0008]
vars = {c[t, x], t, {x}};
RegularizedDeltaPoint[g_, X_List, Xs_List] :=
Piecewise[{{Times @@ Thread[1/(4 g) (1 + Cos[\[Pi]/(2 g) (X - Xs)])],
And @@ Thread[RealAbs[X - Xs] <= 2 g]}, {0, True}}]
Subscript[h, mesh] = Sqrt[Min[mesh["MeshElementMeasure"]]];
Subscript[gamma, reg] = Subscript[h, mesh]/2;
temp = RegularizedDeltaPoint[Subscript[gamma, reg], {x},
includePoints[]];
parameters = {kappa -> {{910}}, v1 -> 162,
gamma -> Subscript[gamma, reg], Qp -> 1.5};
pde = {Derivative[1, 0][c][t, x] +
Inactive[
c[t, x], {x}], {x}] + {v1}.Inactive[Grad][c[t, x], {x}] +
Qp*RegularizedDeltaPoint[gamma, {x}, {10}] == 0,
c[0, x] == 1} /. parameters;

tEnd = 2;
cfun = NDSolveValue[{pde, DirichletCondition[c[t, x] == 5, x == 0]}, c, {t, 0, tEnd}, {x} \[Element] mesh];
Manipulate[
Plot[cfun[t, x], {x} \[Element] region,
PlotRange -> {{0, 100}, {0, 5}}], {t, 0, tEnd}] The above figure displays solution vs x-position. Instead, I want to plot solution curves observed at x-positions in includePoints = {{10}, {20}, {30}, {40}, {50}}; as a function of time. as a function of time.

Suggestions will be really appreciated.

EDIT:

1. May I also know how to save the solution in cfun at includePoints = {{10}, {20}, {30}, {40}, {50}} over the integration time span to a text file?

2. Using the following command

With[{i = Flatten[includePoints]}, Plot[Evaluate[cfun[t, #] & /@ i], {t, 0, tEnd}, PlotRange -> {{0, tEnd}, {0, 5}}, PlotLegends -> (StringTemplate["c(t,)"] /@ i)]]

I could generate, But I am not sure why c(t,0), the first point which is the left boundary isn't at 5 (Dirichlet conditioned defined) at t=0. Could someone please have a look? Also, the legend labels aren't formatted correctly String Template appears as text in the legend).

• Try including the Dirichlet condition in a list with the PDE, like this cfun = NDSolveValue[{pde, DirichletCondition[c[t, x] == 5, x == 0]}, c, {t, 0, tEnd}, {x} \[Element] mesh]; Jun 20 at 8:11
• @LouisB Thank you, I tried the above but it is still not clear to me how the plot can be generated. Could you please explain a bit more? Or, may I know how to save the solution in cfun at includePoints = {{10}, {20}, {30}, {40}, {50}} over the integration time span to a text file? Jun 20 at 10:40
• This question is a continuation of a question here Jun 21 at 4:31

# Update to fix boundary and initial condition inconsistency

Your "MaxCellMeasure" is much more refined than it needs to be. Furthermore, you will note that the initial condition and the DirichletCondition are inconsistent. This causes the DirichletCondition to be ignored. One way to remove this inconsistency is to rapidly ramp up the DirichletCondition from the initial condition to its desired value.

Needs["NDSolveFEM"]
region = Line[{{0}, {100}}];
includePoints = {{10}, {20}, {30}, {40}, {50}};
mesh = ToElementMesh[region, "IncludePoints" -> includePoints,
"MaxCellMeasure" -> 1];
vars = {c[t, x], t, {x}};
RegularizedDeltaPoint[g_, X_List, Xs_List] :=
Piecewise[{{Times @@ Thread[1/(4 g) (1 + Cos[π/(2 g) (X - Xs)])],
And @@ Thread[RealAbs[X - Xs] <= 2 g]}, {0, True}}]
Subscript[h, mesh] = Sqrt[Min[mesh["MeshElementMeasure"]]];
Subscript[gamma, reg] = Subscript[h, mesh]/2;
temp = RegularizedDeltaPoint[Subscript[gamma, reg], {x},
includePoints[]];
parameters = {kappa -> {{910}}, v1 -> 162,
gamma -> Subscript[gamma, reg], Qp -> 1.5};
pde = {Derivative[1, 0][c][t, x] +
Inactive[Div][(-kappa) .
Inactive[Grad][c[t, x], {x}], {x}] + {v1} .
Qp*RegularizedDeltaPoint[gamma, {x}, {10}] == 0, c[0, x] == 1} /.
parameters;

tEnd = 2;
cfun = NDSolveValue[
pde~Join~{DirichletCondition[c[t, x] == 4 (1 - Exp[-1000 t]) + 1,
x == 0]}, c, {t, 0, tEnd}, {x} ∈ mesh];
With[{i = Flatten[{0}~Join~includePoints]},
Plot[Evaluate[cfun[t, #] & /@ i], {t, 0, tEnd},
PlotRange -> {{0, tEnd}, {0, 5.1}},
PlotLegends -> (StringTemplate["c(t,)"] /@ i)]]
ListPlot[Table[cfun[t, x], {x, 0, 100, 10}, {t, 0, tEnd, tEnd/100}],
Joined -> True, DataRange -> {0, tEnd}] Here is one way you could plot the solutions:

With[{i = Flatten[includePoints]},
Plot[Evaluate[cfun[t, #] & /@ i], {t, 0, tEnd},
PlotRange -> {{0, tEnd}, {0, 5}},
PlotLegends -> (StringTemplate["c(t,)"] /@ i)]
] • Thanks so much, works great. Sorry for another question. Could you please have a look at my edit? Jun 20 at 12:44
• @Natasha If I quit the kernel and copy the code from your post, I get the image shown in my answer. Unless you changed the definition of includePoints, 0 should not be there. What \$Version of Mathematica are you using? Jun 20 at 13:04
• Thanks a ton for the answer :) . I'm using version 12.0 Jun 20 at 17:11
• Could you please have a look at this post? Jun 21 at 16:18
• Is there an option to specify MinCellMeasure ? Could you please have a look at this post? Jun 28 at 1:45