# NDSolve Stiff PDE

I'm new to Mathematica and have a hard time solving a stiff PDE using NDSolve after searching around for different stiff system methods. Below is the PDE and so far I have been getting the errors of NDSolve::ndsz: At t == xxx, step size is effectively zero; singularity or stiff system suspected. It's been quite frustrating and any feedback is deeply appreciated.

(*parameter setup*)
pm = 0.64;
p0 = 0.5;
tmax = (pm - p0)/pm;
Tanhm1[t_, t0_] := Tanh[500*(t - t0)] + 1;
T[t_] := 293 + 20*Tanhm1[t, 0.1*tmax];
R = 225*^-9; k = 1.381*^-23;
mu[t_] := (2.414*^-5)*10^(247.8/(T[t] - 140));
g = 9.8; rp = 1050; rl = 1000;
D0[t_] := k*T[t]/(6*Pi*mu[t]*R);
U0[t_] := 2*R^2*g*(rp - rl)/(9*mu[t]);
H0 = 28*^-6; pinf = 0;
A1 = 8.07131; B1 = 1730.63; C1 = 233.426;
pint[t_] := 10^(A1 - B1/(C1 + T[t]));
kg = 20*^-16;
delta[t_] := kg*H0/(rl *D0[t]);
E0[t_] := delta[t]*(pint[t] - pinf);
Pe[t_] := E0[t]*H0/D0[t];
Ns[t_] := U0[t]/E0[t];
K[p_] := (1 - p)^6.55;
Z[p_] := 1/(pm - p);
DP [p_] := D[p*Z[p], p];

(*PDE*)
pde = {D[p[t, x], t] + (x/(1 - t))*D[p[t, x], x] ==
Ns[t]/(1 - t)*D[K[p[t, x]]*p[t, x], x] +
1/(Pe[t]*(1 - t)^2)*D[K[p[t, x]]*DP[p[t, x]]*D[p[t, x], x], x]};
bc = {Ns[t]*Pe[t]*p[t, 0]*(1 - t)*(1 - Exp[-1000*t]) +
DP[p[t, 0]]*Derivative[0, 1][p][t, 0] == 0,
Ns[t]*Pe[t]*K[p[t, 1]]*p[t, 1]*(1 - t)*(1 - Exp[-1000*t]) +
K[p[t, 1]]*DP[p[t, 1]]*Derivative[0, 1][p][t, 1] ==
Pe[t]*(1 - t)*p[t, 1]*(1 - Exp[-1000*t]),
p[0, x] == p0};
sol = NDSolve[Flatten[{pde, bc}], p, {t, 0, tmax}, {x, 0, 1},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 41, "MaxPoints" -> 101,
"DifferenceOrder" -> "Pseudospectral"}}]
Plot3D[Evaluate[p[t, x] /. sol], {t, 0, tmax}, {x, 0, 1}, PlotPoints -> 100]
$$$$

• Removing "DifferenceOrder" -> "Pseudospectral" seems to help. Also, add the option PlotRange -> All in Plot3D to see the entire solution. Commented Mar 14, 2020 at 14:38
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) 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, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Commented Mar 14, 2020 at 14:40

Removing all the options to NDSolve, it works fine in version 8.0. Simpliy do

sol = NDSolve[Flatten[{pde, bc}], p, {t, 0, tmax}, {x, 0, 1}]

Plot3D[Evaluate[p[t, x] /. sol], {t, 0, tmax}, {x, 0, 1},
PlotPoints -> 100, PlotRange -> All]


Can you use the StiffnessSwithching method?

 sol = NDSolve[Flatten[{pde, bc}], p, {t, 0, tmax}, {x, 0, 1},
Method -> {"StiffnessSwitching"}]

Plot3D[Evaluate[p[t, x] /. sol], {t, 0, tmax}, {x, 0, 1},
PlotPoints -> 100, PlotRange->All]
`
• Sorry, @maeinss, i edited your code accidentally when trying to edit mine. Don't know how to undo. Commented Apr 15, 2020 at 6:43