# Where is the numerical solving breaking down?

I am working with a set of three coupled reaction-diffusion PDEs, and for some parameter values I am getting some not so great solutions. I have been searching documentation and tutorials, and I have been trying to implement various methods and change various settings, but I feel like I am just blindly trying things without fully understanding what is going on. The equations might look intimidating, but I have supplied my code as well to help make everything reproducible.

As a summary, I was asking about why a bad solution is obtained for a certain parameter set. As I went through testing my problem to make this post I came across a second issue where the solution behavior also depended on the end time of integration, which seems odd to me. Why am I running into these issues?

# Equations

The set of equations I am working with is as follows: $$\frac{\partial x}{\partial t}=\left(1+\frac{\gamma_{pL}\cdot y_T}{(1+\gamma_P\cdot x)^2}\right)^{-1}\cdot\left[\gamma_{D_P}\cdot\nabla^2x+\gamma_a\cdot x_I+\gamma_c\cdot\frac{\gamma_{pL}}{\gamma_P}\cdot\left(\frac{\gamma_P\cdot x}{1+\gamma_P\cdot x}\right)^2\cdot y_T\right]$$

$$\frac{\partial y_T}{\partial t}=-\gamma_c\cdot\frac{\gamma_P\cdot x}{1+\gamma_P\cdot x}\cdot y_T$$

$$\frac{\partial z}{\partial t}=\left(1+\frac{\gamma_R}{(1+\gamma_L\cdot z)^2}\right)^{-1}\cdot\left[\gamma_{D_L}\cdot\nabla^2z+\gamma_c\cdot\frac{\gamma_P\cdot x}{1+\gamma_P\cdot x}\cdot y_T\right]$$

where $$x$$, $$y_T$$, and $$z$$ are functions of $$r$$ (polar coordinate) and $$t$$ (time), $$x_I$$ is a source term given by $$x_I=\frac{1}{2\cdot\gamma_{D_P}\cdot t+1}\cdot\exp\left(-\frac{r^2}{2\cdot(2\cdot\gamma_{D_P}\cdot t+1)}-\gamma_a\cdot t\right)$$ and all $$\gamma_i$$ terms are parameters.

Initial conditions and boundary conditions are as follows: $$x(r,0)=z(r,0)=0,\ y_T(r,0)=1$$ $$\left.\frac{\partial x}{\partial t}\right|_{r=0}=\left.\frac{\partial z}{\partial t}\right|_{r=0}=0$$ $$x(\infty,0)=z(\infty,0)=0$$

In words, we start off with no $$x$$ and $$z$$, and $$y_T$$ starts at a maximum value of $$1$$ across all space. The no-flux boundary condition at the origin is because the system has radial symmetry, and we want $$x$$ and $$z$$ to go to $$0$$ at $$\infty$$.

# Code

To tweak the above equations to work nicely in Mathematica (I think?), we need to use a "large" value of $$r$$ in place of $$\infty$$ and we need to use a point close to $$r=0$$ rather than $$r=0$$ itself (due to the $$1/r$$ term in the Laplacian in polar coordinates).

dr = .001; (*"small" r since we cannot use the origin in polar coordinates*)
ρFar = 100.; (*r that is "infinity"*)
τmax = 6.38;(*time to solve out to*)


The equations above:

xI = 1/(2*γDP*t + 1)*Exp[-(r-dr)^2/(2*(2*γDP*t + 1)) - γa*t];

dx = D[x[r, t],t] == (1 + (γpL*yT[r, t])/(1 + γP*x[r,t])^2)^-1*(γDP*Laplacian[x[r, t], {r, θ}, "Polar"] + γa*xI + γc*γpL/γP*((γP*x[r, t])/(1 + γP*x[r, t]))^2*yT[r, t]);

dyT = D[yT[r, t], t] == (-γc*γP*x[r, t]*yT[r, t])/(1 + γP*x[r, t]);

dz = D[z[r, t],t] == (1 + γR/(1 + γL*z[r, t])^2)^-1*(γDL*Laplacian[z[r, t], {r, θ},"Polar"] + (γc*γP*x[r, t]*yT[r, t])/(1 + γP*x[r, t]));


Intial/Boundary Conditions:

initx = x[r, 0] == 0.;

bc1x = (D[x[r, t], r] == 0) /. r -> dr;

bc2x = x[ρFar, t] == 0.;

inityT = yT[r, 0] == 1.;

initz = z[r, 0] == 0;

bc1z = (D[z[r, t], r] == 0) /. r -> dr;

bc2z = z[ρFar, t] == 0.;


Putting it all together:

deqns = {dx, dyT, dz, initx, bc1x, bc2x, inityT, initz, bc1z, bc2z};


### Problem

My problem arises for certain parameter values (not all). I have also found that using

MaxStepSize -> .1


for NDSolveValue seems to help with some instances of my problem, but not all. An example of a "good solution" can be obtained with the following parameters

paramsGood = {γDP -> 0.001222, γDL -> 122.2, γa -> 1., γc -> 100., γP -> 0.01, γpL -> 0.01, γL -> 100., γR -> 0.01};

zSoln = NDSolveValue[deqns /. paramsGood, z, {r, dr, ρFar}, {t,0., τmax}, MaxStepSize -> .1];

Plot[zSoln[r, 2.6], {r, dr, 15}, PlotRange -> {-.002, .002}]


However, the following bad parameter set

paramsBad = {γDP -> 0.1222, γDL -> 122.2, γa -> 100., γc -> 1., γP -> 0.01, γpL -> 100., γL -> 1000., γR -> 100.};


gives an issue in the solution for $$z$$ (I have picked the time $$t=2.6$$ to show the bad solution, but the solution doesn't seem good from the start).

zSoln = NDSolveValue[deqns /. paramsBad, z, {r, dr, ρFar}, {t,0., τmax}, MaxStepSize -> .1];
Plot[zSoln[r, 2.6], {r, dr, 15}, PlotRange -> {-.002, .002}]


None of the dependent variables should be negative.

I have tried modifying the MaxStepSize, StartingStepSize, and specifying a spatial discretization as suggested in the documentation:

Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 1000}}}


but I cannot seem to fix this issue. Sometimes these things fix the problem, but then for other parameter values the issues come right back (for example, I think making the MinPoints option above larger fixes for this parameter set, but not for all of the. Additionally the solution ends up taking a long time). Why is this happening for the numerical solver, and what is the most efficient way to determine a well-behaved solution Note that I need to be able to solve these equations fairly quickly.

### Second Problem

If I just change the end time to $$t=6.5$$ the solution ends up behaving differently. There is a very small response in $$z$$, but it at least looks well-behaved?

zSoln = NDSolveValue[deqns /. paramsBad, z, {r, dr, ρFar}, {t,0., 6.5}, MaxStepSize -> .1];
Plot[zSoln[r, 2.6], {r, dr, 15}, PlotRange -> {-.002, .002}]


Also, this issue does not occur for the "good parameter set"; the solution does not show any notable difference in choice of the final integration time.

So this really makes me question how do I make sure I am getting a legitimate solution here? Which settings should I employ to be confident in my solutions for a wide range of parameters?

• MaxStepSize -> .1 is not close to providing adequate resolution near t = 0 and r = 0 for paramsBad. In fact it really is not adequate for paramsGood. However, just decreasing MaxStepSize to .001 leads to a very slow computation. I am exploring various options. I should add that you really need to use Plot3D to see whether your answers are reasonable. – bbgodfrey Oct 18 '19 at 4:36
• Also, running to t = 6.5 does not solve the problem either. It just moves it to a different time slice. "Why" remains a good question, though. – bbgodfrey Oct 18 '19 at 4:44
• @bbgodfrey Thanks for looking into it. I tend to use the Manipulate function to move through time rather than Plot3D. But rest assured I'm not just looking at single times. – Aaron Stevens Oct 18 '19 at 4:52
• Look very near t = 0, r = 0 to see a negative spike even for tmax = 6.5 This differs from tmax = 6.38, because the temporal zoning changes. Use zSoln["Coordinates"] to see the zoning.. I shall look at this some more tomorrow. – bbgodfrey Oct 18 '19 at 5:07

Some experimentation shows that the NDSolve problem is associated exclusively with z with no feedback to x or yT, because dz and dyT are independent of z. It is, therefore, convenient first to compute x and yT, and only then to compute z. For paramsGood,

deqns = {dx, dyT, initx, bc1x, bc2x, inityT};
deqnz = {dz, initz, bc1z, bc2z};

{xSoln, yTSoln} = NDSolveValue[deqns /. paramsGood, {x, yT}, {r, dr, ρFar},
{t, 0., τmax}, MaxStepSize -> .1];
zSoln = NDSolveValue[deqnz /. paramsGood /. {x[r, t] -> xSoln[r, t], yT[r, t] -> yTSoln[r, t]},
z, {r, dr, ρFar}, {t, 0., τmax}, MaxStepSize -> .1];
Plot3D[xSoln[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All, AxesLabel -> {r, t, x},
PlotPoints -> 50, LabelStyle -> {Bold, Black, 15}, ImageSize -> Large,
ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[yTSoln[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All, AxesLabel -> {r, t, y},
PlotPoints -> 50, LabelStyle -> {Bold, Black, 15}, ImageSize -> Large,
ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[zSoln[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All, AxesLabel -> {r, t, z},
PlotPoints -> 50, LabelStyle -> {Bold, Black, 15}, ImageSize -> Large,
ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]


Visibly, the three computed veriables are well behaved everywhere. A slice through z at t = 2.6 reproduces the first plot in the question, as expected. The same code applied to paramsBad (but with PlotPoints -> 250) does not do so well.

Note the spikes near the origin in the plot of z. (The z computation in this case also produces a few underflow warnings, which I believe are insignificant.) Several constant-r slices through the z plot for very small r yield.

Plot[Table[zSoln[0.01 n, t], {n, 1, 9, 2}], {t, 0, .04}, PlotRange -> All, AxesLabel -> {t, z}]


On the other hand, a plot corresponding to the second plot in the question looks fine, although it may not be.

Plot[{zSoln[r, 2.6]}, {r, dr, 15}, PlotRange -> All, AxesLabel -> {r, z}]


In any case, the solution for z obtained here for paramsBad is much better than that in the question, probably because separating the computations of z allows NDSolve to optimize its automatic time-step determination for z alone.

Finally, let us repeat the computation but with MaxStepSize -> .01. Plots for x and yT appear identical to those just given, with the z plot now seemingly without error.

But, it isn't.

Plot[Table[zSoln[0.01 n, t], {n, 1, 9, 2}], {t, 0, .04}, PlotRange -> All, AxesLabel -> {t, z}]


Better, but not perfect. Note that the slice for t = 2.6 changes a bit in this case.

To answer the specific questions posed by the OP, first the errors obtained for paramsBad are due to lack of resolution, largely corrected by reducing MaxStepSize to 0.01. Still smaller MaxStepSize undoubtedly would produce a still more accurate result, but the computation becomes quite slow for 0.001. I also experimented some with very small StartingStepSize, but obtained little improvement. Second, the apparent improvement from the second to the third plots in the question are due to different time-steps used by NDSolve in the two cases, but plotting z in 3D shows that significant errors persist in the computation. Therefore, little credence should be given to the third plot in the question.

• Thank you for a more thorough look into this. So then, what is my best course of action here? It looks like Solving equations for $x$ and $y_T$ individually first might help a bit, and then reducing the MaxStepSize I suppose? – Aaron Stevens Oct 18 '19 at 18:12
• Yes, that is correct. – bbgodfrey Oct 18 '19 at 18:32
• Do you think it would be worth it to put a bounty and see if other people have suggestions? Or do you think what you have suggested is as far as this can go? I just don't know what else to try. – Aaron Stevens Oct 19 '19 at 0:31
• Who knows? @MichaelE2 might have some ideas. – bbgodfrey Oct 19 '19 at 1:44

For this problem, you must specify a solution method. Since the system of equations is nonlinear and the equation for y does not contain derivatives with respect to r, we must choose "MethodOfLines". It solves all problems right away. The code contains Quiet, since function xI is used outside the machine number, for example, 1/2.71828^726.112is calculated.

dr = .001;(*"small" r since we cannot use the origin in polar \
coordinates*)ρFar = 100.;(*r that is "infinity"*)τmax = \
6.38;(*time to solve out to*)
xI = 1/(2*γDP*t + 1)*
Exp[-(r - dr)^2/(2*(2*γDP*t + 1)) - γa*t];

dx = D[x[r, t],
t] == (1 + (γpL*
yT[r, t])/(1 + γP*x[r, t])^2)^-1*(γDP*
Laplacian[x[r, t], {r, θ}, "Polar"] + γa*
xI + γc*γpL/γP*((γP*
x[r, t])/(1 + γP*x[r, t]))^2*yT[r, t]);

dyT = D[yT[r, t],
t] == (-γc*γP*x[r, t]*
yT[r, t])/(1 + γP*x[r, t]);

dz = D[z[r, t],
t] == (1 + γR/(1 + γL*z[r, t])^2)^-1*(γDL*
Laplacian[z[r, t], {r, θ},
"Polar"] + (γc*γP*x[r, t]*
yT[r, t])/(1 + γP*x[r, t]));
initx = x[r, 0] == 0.;

bc1x = (D[x[r, t], r] == 0) /. r -> dr;

bc2x = x[ρFar, t] == 0.;

inityT = yT[r, 0] == 1.;

initz = z[r, 0] == 0;

bc1z = (D[z[r, t], r] == 0) /. r -> dr;

bc2z = z[ρFar, t] == 0.;
deqns = {dx, dyT, dz, initx, bc1x, bc2x, inityT, initz, bc1z, bc2z};
paramsGood = {γDP -> 0.001222, γDL ->
122.2, γa -> 1., γc -> 100., γP ->
0.01, γpL -> 0.01, γL -> 100., γR -> 0.01};

zSoln = NDSolveValue[deqns /. paramsGood,
z, {r, dr, ρFar}, {t, 0., τmax},
Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> "Residual",
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 137, "MaxPoints" -> 137,
"DifferenceOrder" -> 2}}}] // Quiet;

paramsBad = {γDP -> 0.1222, γDL -> 122.2, γa ->
100., γc -> 1., γP -> 0.01, γpL ->
100., γL -> 1000., γR -> 100.};
z, {r, dr, ρFar}, {t, 0., τmax},
Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> "Residual",
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 137, "MaxPoints" -> 137,
"DifferenceOrder" -> 2}}}] // Quiet;

z, {r, dr, ρFar}, {t, 0., 6.5},
Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> "Residual",
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 137, "MaxPoints" -> 137,
"DifferenceOrder" -> 2}}}] // Quiet;

{Plot3D[zSoln[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All,
PlotLabel -> "paramsGood", AxesLabel -> Automatic, Mesh -> None,
ColorFunction -> "Rainbow"],
Plot3D[zSoln1[r, t], {r, dr, 15}, {t, 0, τmax},
PlotRange -> All, PlotLabel -> "paramsBad", AxesLabel -> Automatic,
ColorFunction -> "Rainbow"],
Plot3D[zSoln2[r, t], {r, dr, 15}, {t, 0, 6.5}, PlotRange -> All,
PlotLabel -> "paramsBad", AxesLabel -> Automatic,
ColorFunction -> "Rainbow"]}


Consider the following example

paramsGood = {γDP -> 0.01222, γDL -> 122.2, γa ->
0.01, γc -> 100., γP -> 1., γpL ->
1., γL -> 1000., γR -> 100.}


In this case, for "MaxPoints" -> 137, the solution goes into region z<0 (Fig. 2, the upper row on the left). We increase the number of points to 1037 and to 2037 (the top row is the center and to the right). The solution converges. However, with an increase of "DifferenceOrder" ->4 , the solution looks different (bottom row).This is due to artificial viscosity, which decreases with increasing "DifferenceOrder". I can recommend a combination of options "MinPoints" -> 1037, "MaxPoints" -> 1037, "DifferenceOrder" -> 2 for such cases (with the number of points you can experiment).

• Thank you for this. I had been trying to specify the method, but I guess I was not being specific enough in setting aspects of the method. Which method settings here are the best to tweak if I am still getting some poor solutions? It looks like fewer of my parameter tests are producing bad solutions with your method specification, but I am still getting some bad solutions for other parameter values. – Aaron Stevens Oct 19 '19 at 18:43
• And I will add that some of the break downs do occur at $r$ values farther from $r=0$ (like from $r\approx2$ to $r\approx8$) – Aaron Stevens Oct 19 '19 at 18:52
• @AaronStevens Can you write something specific about the parameters. – Alex Trounev Oct 19 '19 at 18:53
• {\[Gamma]DP -> 0.01222, \[Gamma]DL -> 122.2, \[Gamma]a -> 0.01, \[Gamma]c -> 100., \[Gamma]P -> 1., \[Gamma]pL -> 1., \[Gamma]L -> 1000., \[Gamma]R -> 100.} I think I have been doing my tests with the same code as in my question. Let me know if everything actually works fine for those parameters. I have been doing many various tests since posting :) Also I removed the MaxStep specification that I had in my post, since your answer does not have it. – Aaron Stevens Oct 19 '19 at 18:55
• @AaronStevens Sorry, fixed it on "DifferenceOrder"->2 – Alex Trounev Oct 19 '19 at 22:52

# Fully NDSolve-based Solution

Adjustion for spatial step size together with temporal step size helps. I've used parameters mentioned in the comment for testing:

mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}

paramsBad = {γDP -> 0.01222, γDL -> 122.2, γa -> 0.01, γc ->
100., γP -> 1., γpL -> 1., γL -> 1000., γR -> 100.};

difforder = 4; points = 200;
test = NDSolveValue[
deqns /. paramsBad, {x, yT, z}, {r, dr, ρFar}, {t, 0, τmax},
Method -> mol[points, difforder], MaxStepSize -> {Automatic, 0.01}];
Plot3D[test[[-1]][r, t], {r, eps, 15}, {t, 0, τmax}, PlotRange -> All]


# Partly NDSolve-based Solution

Another possible solution is discretizing the system in $$r$$ direction all by ourselves, in this case NDSolve seems to do a better job on choosing temporal step size so we don't need to adjust MaxStepSize option. Also, the singularity at $$r=0$$ is removed i.e. we can use dr = 0 in this method.

Notice pdetoode is used for discretization in the following program.

Clear[dr]
ρFar = 100.;
τmax = 6.38;
xI = 1/(2 γDP t + 1) Exp[-(r - dr)^2/(2 (2 γDP t + 1)) - γa t];
dx = D[x[r, t], t] == (1 + (γpL yT[r, t])/(1 + γP x[r, t])^2)^-1 (γDP
Laplacian[x[r, t], {r, θ}, "Polar"] + γa
xI + γc γpL/γP ((γP x[r, t])/(1 + γP x[r, t]))^2 yT[r, t]);
dyT = D[yT[r, t], t] == (-γc γP x[r, t] yT[r, t])/(1 + γP x[r, t]);
dz = D[z[r, t], t] == (1 + γR/(1 + γL z[r, t])^2)^-1 (γDL
Laplacian[z[r, t], {r, θ}, "Polar"] + (γc γP x[r, t] yT[r, t])/(1 + γP x[r, t]));
initx = x[r, 0] == 0.;
bc1x = (D[x[r, t], r] == 0) /. r -> dr;
bc2x = x[ρFar, t] == 0.;
inityT = yT[r, 0] == 1.;
initz = z[r, 0] == 0;
bc1z = (D[z[r, t], r] == 0) /. r -> dr;
bc2z = z[ρFar, t] == 0.;
removesingularity = eq \[Function] r # & /@ Expand[eq] // Expand;
{newdx, newdz} = removesingularity /@ {dx, dz};
domain = {dr, ρFar};
points = 100; difforder = 2; grid = Array[# &, points, domain];
ptoofunc = pdetoode[{x, yT, z}[r, t], t, grid, difforder];
del = #[[2 ;; -2]] &;
ode = {del /@ ptoofunc@{newdx, newdz}, ptoofunc@dyT};
odeic = ptoofunc@{initx, inityT, initz};
odebc = With[{sf = 1}, ptoofunc@Map[sf # + D[#, t] &, {bc1x, bc2x, bc1z, bc2z}, {2}]];
var = Outer[#[#2] &, {x, yT, z}, grid];
paramsBad = {γDP -> 0.01222, γDL -> 122.2, γa -> 0.01, γc ->
100., γP -> 1., γpL -> 1., γL -> 1000., γR -> 100.};
eps = 0;
sollst = NDSolveValue[{ode, odeic, odebc} /. paramsBad /. dr -> eps,
var /. dr -> eps, {t, 0, τmax}];
{solx, soly, solz} = rebuild[#, grid /. dr -> eps, 2] & /@ sollst;
Plot3D[solz[r, t], {r, 0, 15}, {t, 0, τmax}, PlotRange -> All]


## Remark

1. As one can see, only 100 grid points are used. Increasing points won't cause significant difference.

2. Though difforder = 2 is chosen, difforder = 4 can be used. You may need to add Method -> {"EquationSimplification" -> Solve} option to NDSolve though.

3. Using dr = 0.001 won't cause any significant difference.

I think a key issue is the inability to resolve your source term, $$x_I$$, relative to your domain. Three parameter sets have been given through the course of the discussion. If we assume that we will use the same spatial and temporal domain for each parameter set, then a Plot3D of the source term reveals the relative feature size of the pulse.

Plot3D[xI /. paramsGood, {t, 0, τmax}, {r, dr, ρFar},
PlotRange -> All, PlotLabel -> "Source Term paramsGood",
AxesLabel -> Automatic, Mesh -> None, ColorFunction -> "DarkBands",
PlotPoints -> 50, MaxRecursion -> 4]
Plot3D[xI /. paramsBad, {t, 0, τmax}, {r, dr, ρFar},
PlotRange -> All, PlotLabel -> "Source Term paramsBad",
AxesLabel -> Automatic, Mesh -> None, ColorFunction -> "DarkBands",
PlotPoints -> 50, MaxRecursion -> 4]
Plot3D[xI /. paramsBad2, {t, 0, τmax}, {r, dr, ρFar},
PlotRange -> All, PlotLabel -> "Source Term paramsBad2",
AxesLabel -> Automatic, Mesh -> None, ColorFunction -> "DarkBands",
PlotPoints -> 50, MaxRecursion -> 4]


As you can see, the first bad parameter set operates on a tiny portion of the domain. When the ratio of largest relevant feature to the smallest relevant feature that you want to capture is large, discretization can be difficult. In this case, we can use a helper function to create a non-uniform mesh that is very fine near $$r=0$$ and expands 100,000 times larger at the end of the domain with 200 grid points.

meshGrowth[x0_, xf_, n_, ratio_] := Module[{k, fac, delta},
k = Log[ratio]/(n - 1);
fac = Exp[k];
delta = (xf - x0)/Sum[fac^(i - 1), {i, 1, n - 1}];
N[{x0}~Join~(x0 +
delta Rest@
FoldList[(#1 + #2) &, 0,
PowerRange[fac^0, fac^(n - 3), fac]])~Join~{xf}]
]
grid = meshGrowth[0.001, 100, 200, 100000];
ListLogPlot@grid


We can now apply this non-uniform grid to NDSolve. We also note that the source pulse is very short lived for the first parameters set. We will set the StartingStepSize and MaxStepFraction to small values to resolve the time domain. If we apply this to @bbgodfrey's approach, we obtain the following:

paramsGood = {γDP -> 0.001222, γDL ->
122.2, γa -> 1., γc -> 100., γP ->
0.01, γpL -> 0.01, γL -> 100., γR -> 0.01};
paramsBad = {γDP -> 0.1222, γDL -> 122.2, γa ->
100., γc -> 1., γP -> 0.01, γpL ->
100., γL -> 1000., γR -> 100.};
paramsBad2 = {γDP -> 0.01222, γDL -> 122.2, γa ->
0.01, γc -> 100, γP -> 1, γpL ->
1., γL -> 1000., γR -> 100.};
dr = .001;(*"small" r since we cannot use the origin in polar \
coordinates*)
ρFar = 100.;(*r that is "infinity"*)
τmax = 6.38;(*time to solve out to*)
opts = (Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"Coordinates" -> {grid}}});
{xSolnG, yTSolnG} =
NDSolveValue[
deqns /. paramsGood, {x, yT}, {r, dr, ρFar}, {t,
0., τmax}, opts, StartingStepSize -> 0.0001,
MaxStepFraction -> 1/400];
zSolnG = NDSolveValue[
deqnz /. paramsGood /. {x[r, t] -> xSolnG[r, t],
yT[r, t] -> yTSolnG[r, t]},
z, {r, dr, ρFar}, {t, 0., τmax}, opts,
StartingStepSize -> 0.0001, MaxStepFraction -> 1/400];
Plot3D[xSolnG[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All,
AxesLabel -> {r, t, x}, PlotPoints -> 50,
LabelStyle -> {Bold, Black, 15}, ImageSize -> Large,
ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[yTSolnG[r, t], {r, dr, 15}, {t, 0, τmax},
PlotRange -> All, AxesLabel -> {r, t, y}, PlotPoints -> 50,
LabelStyle -> {Bold, Black, 15}, ImageSize -> Large,
ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[zSolnG[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All,
AxesLabel -> {r, t, z}, PlotPoints -> 50,
LabelStyle -> {Bold, Black, 15}, ImageSize -> Large,
ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]

{xSolnB, yTSolnB} =
NDSolveValue[
deqns /. paramsBad, {x, yT}, {r, dr, ρFar}, {t,
0., τmax}, opts, StartingStepSize -> 0.0001,
MaxStepFraction -> 1/400];
zSolnB = NDSolveValue[
deqnz /. paramsBad /. {x[r, t] -> xSolnB[r, t],
yT[r, t] -> yTSolnB[r, t]},
z, {r, dr, ρFar}, {t, 0., τmax}, opts,
StartingStepSize -> 0.0001, MaxStepFraction -> 1/400];
Plot3D[xSolnB[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All,
AxesLabel -> {r, t, x}, PlotPoints -> 50,
LabelStyle -> {Bold, Black, 15}, ImageSize -> Large,
ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[yTSolnB[r, t], {r, dr, 15}, {t, 0, τmax},
PlotRange -> All, AxesLabel -> {r, t, y}, PlotPoints -> 50,
LabelStyle -> {Bold, Black, 15}, ImageSize -> Large,
ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[zSolnB[r, t], {r, dr, 15}, {t, 0, τmax}, PlotRange -> All,
AxesLabel -> {r, t, z}, PlotPoints -> 50,
LabelStyle -> {Bold, Black, 15}, ImageSize -> Large,
ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]

{xSolnB2, yTSolnB2} =
NDSolveValue[
deqns /. paramsBad2, {x, yT}, {r, dr, ρFar}, {t,
0., τmax}, opts, StartingStepSize -> 0.0001,
MaxStepFraction -> 1/400];
zSolnB2 =
NDSolveValue[
deqnz /. paramsBad2 /. {x[r, t] -> xSolnB2[r, t],
yT[r, t] -> yTSolnB2[r, t]},
z, {r, dr, ρFar}, {t, 0., τmax}, opts,
StartingStepSize -> 0.0001, MaxStepFraction -> 1/400];
Plot3D[xSolnB2[r, t], {r, dr, 15}, {t, 0, τmax},
PlotRange -> All, AxesLabel -> {r, t, x}, PlotPoints -> 50,
LabelStyle -> {Bold, Black, 15}, ImageSize -> Large,
ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[yTSolnB2[r, t], {r, dr, 15}, {t, 0, τmax},
PlotRange -> All, AxesLabel -> {r, t, y}, PlotPoints -> 50,
LabelStyle -> {Bold, Black, 15}, ImageSize -> Large,
ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]
Plot3D[zSolnB2[r, t], {r, dr, 15}, {t, 0, τmax},
PlotRange -> All, AxesLabel -> {r, t, z}, PlotPoints -> 50,
LabelStyle -> {Bold, Black, 15}, ImageSize -> Large,
ViewPoint -> {-1, -3, 5/4}, ViewVertical -> {1/8, 1/3, 1}]

Plot[Table[zSolnG[0.01 n, t], {n, 1, 9, 2}], {t, 0, .04},
PlotRange -> All, AxesLabel -> {t, z}, PlotLabel -> "Good"]
Plot[Table[zSolnB[0.01 n, t], {n, 1, 9, 2}], {t, 0, .04},
PlotRange -> All, AxesLabel -> {t, z}, PlotLabel -> "Bad"]
Plot[Table[zSolnB2[0.01 n, t], {n, 1, 9, 2}], {t, 0, .04},
PlotRange -> All, AxesLabel -> {t, z}, PlotLabel -> "Bad2"]


We can also apply this non-uniform mesh to you original set of coupled equations for the first bad (i.e., sharpest pulse) parameter set and it seems to still hold together.

deqns = {dx, dyT, dz, initx, bc1x, bc2x, inityT, initz, bc1z, bc2z};
zSolnorig =
`