# NDSolve returns solution with single point domain

Thanks to helpful comments from Michael E2 and George2079, I was able to focus in on exactly the source of the issue. With some simplification, I can reduce the problem to:

NDSolve[{Derivative[1][y][t] == (-2 E^(1/5 y[t]) + 1)/(-600 E^(1/5 y[t]) +
20338 + (2/10) t ), y[101690] ==  0}, y, {t, 0, 101690},
MaxSteps -> \[Infinity], WorkingPrecision -> 150]


It quickly throws:

NDSolve::ndcf: Repeated convergence test failure at t == <<185>>; unable to continue. >>


But this problem has an analytical solution:

y[t_] := 1/600 (2 t + C[1] + 3000
ProductLog[C[2],
-((E^(-(t/1500) - C[1]/3000) (101690 + t))/3000)])/.
{C[1] -> 0, C[2] -> -1}


Which you can verify via:

In[1]:= N[y[101690]]
Out[1]= 0.

In[2]:= Simplify[y'[t] -
(-2 E^(1/5 y[t]) + 1)/(-600 E^(1/5 y[t]) + 20338 + (2/10) t)]
Out[2]= 0


This function is pretty well behaved.

I am hoping to understand the related issues with NDSolve, because I'm relying on that engine for more complex problems and this issue is a barrier to its application.

The documentation for the ndcf error says:

If you see this message and the cause is not evident from the nature of the example, please send your example to Technical Support so that it can be investigated.

• I get a step-size-too-small error. You have to increase the precision to around 100 to get effectively nonzero step sizes. But they're still so small (~10^-80 - 10^-30 for PrecisionGoal -> 8) that integration will never finish (effectively). – Michael E2 May 18 '16 at 2:38
• wrapping everything in SetPrecision is unintentionally also converting integers to reals. Now your equations are expressed in terms of variables such as Cee[4.000000000000000000....]. Since now the symbols in the vars list (eg Cee[4]) do not appear in the equations you get a trivial solution. – george2079 May 18 '16 at 20:27
• You can fix that issue by instead doing InitialConditions = {Cee[4][t0] == -270917. ... } /. x_Real :> SetPrecision[x, prec]. (With that fix I get an out of memory error though) – george2079 May 18 '16 at 20:32

The simplified system makes me forget why I suggested the issue concerned working precision. (The edit history shows the OP was much more complicated.) Aside from the fact that @xzczd's answer shows a high WorkingPrecision is not needed, the convergence problem arises from PrecisionGoal and AccuracyGoal being automatically raised (to half WorkingPrecision). That imposes a rather strict requirement with WorkingPrecision -> 150, which apparently NDSolve cannot satisfy. There are two principal reasons for raising WorkingPrecision (WP):

• One cannot achieve the precision and accuracy goals at the current WP. In such a case, raise WP but also manually set PrecisionGoal and AccuracyGoal to their current values (don't raise them).

• One wants to raise PrecisionGoal and AccuracyGoal.

The present case is the first one. We should set PrecisionGoal -> 8 and AccuracyGoal -> 8 (or other desired values). With the high WP, NDSolve can get quite close to the singularity, making quite small steps. I suggest setting MaxSteps to a finite value; 100000 takes a few seconds. You could set MaxSteps -> 10000 (the default) to get the idea.

time = 0; steps = 0;
Dynamic@{time, steps}  (* monitor progress; it's educational *)

steps = 0;
{sol} = NDSolve[{Derivative[1][y][t] == (-2 E^(1/5 y[t]) + 1)/(-600 E^(1/5 y[t]) +
20338 + (2/10) t), y[101690] == 0}, y, {t, 0, 101690},
MaxSteps -> 100000, WorkingPrecision -> 150, PrecisionGoal -> 8,
AccuracyGoal -> 8, StartingStepSize -> 0.001,
StepMonitor :> (time = t; steps++)]

ListLinePlot[N[y /. sol], PlotRange -> All]


Without the 2 options, NDSolve manages to give the correct solution:

nsol = NDSolveValue[{Derivative[1][y][
t] == (-2 E^(1/5 y[t]) + 1)/(-600 E^(1/5 y[t]) + 20338 + (2/10) t), y[101690] == 0},
y, {t, 0, 101690}]


NDSolveValue::ndsz

NDSolve spits out ndsz warning and stops at about 6883.23, but this behavior is reasonable, because when t < 6883.23 there's no real solution:

sol[t_] = 1/
600 (2 t + C[1] +
3000 ProductLog[C[2], -((E^(-(t/1500) - C[1]/3000) (101690 + t))/3000)]) /. {C[
1] -> 0, C[2] -> -1};

{{tL, tR}} = nsol["Domain"];
Plot[{Re@#, Im@#} &@sol@t // Evaluate, {t, 0, tR}, PlotRange -> All,
Epilog -> {Dashed, Line[{{tL, -10}, {tL, 20}}]}]~Show~
Plot[nsol[t], {t, tL, tR}, PlotStyle -> {Thick, Dashed, Red}, PlotRange -> All]