differential equations - Why does this stiff BVP return unevaluated but a similar IVP produce a solution?

In Solving a steady-state viscous Burger's equation with NDSolve, one question involves the following problem:

NDSolve[{u''[x] - 20 u[x]*u'[x] == 0, u[-1] == 1.01, u[1] == -1},
 u[x], {x, -1, 1}]

NDSolve::ndsz: At x == -0.74236, step size is effectively zero; singularity or stiff system suspected.

(* NDSolve[{-20. u[x] Derivative[1][u][x] + (u^\[Prime]\[Prime])[x] == 0,
   u[-1] == 1.01, u[1] == -1}, u[x], {x, -1, 1}] *)

The NDSolve call returns unevaluated. However, the following IVP reports the same problem but produces a partial solution:

NDSolve[{u''[x] - 20 u[x]*u'[x] == 0,
   u[-1] == 1.01, u'[-1] == 0.235320844120522}, 
 u[x], {x, -1, 1}]

NDSolve::ndsz: At x == -0.74236, step size is effectively zero; singularity or stiff system suspected.

What's going on?

Answer

I think the first behavior is arguably a bug and I reported a suggested fix ([CASE:4344046]).

When NDSolve uses the Shooting method, sometimes the initial conditions it tries leads to stiffness. At that point the method fails, and NDSolve returns unevaluated. The only message the user sees the stiffness/singularity message NDSolve::ndsz. However, the real error lies in the initial conditions chosen in by the Shooting method, and that is not reported to the user.

I suggest that the shooting method could check for this, and report that specifying "StartingInitialConditions" might help.

Another possible fix is to tweak the implicit solver used to shoot for the initial conditions to see if we can keep it from overshooting too far. For instance, the first attempt succeeds but poorly. However, we can use it to get good guesses for the initial conditions:

sol = NDSolveValue[
   {-20 u[x] u'[x] + u''[x] == 0, u[-1] == 1.01`, u[1] == -1}, 
   u, {x, -1, 1}, 
   Method -> {"Shooting", 
     "ImplicitSolver" -> {"Newton", 
       "StepControl" -> {"LineSearch", 
         "MaxRelativeStepSize" -> 1/110}}}];

NDSolveValue::berr: The scaled boundary value residual error of 1307.7304404133247` indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found.

(* use sol[-1] and sol[1] for the ICs and recompute *)
sol = NDSolveValue[
   {-20 u[x] u'[x] + u''[x] == 0, u[-1] == 1.01`, u[1] == -1}, 
   u, {x, -1, 1},
   Method -> {"Shooting",
     "StartingInitialConditions" -> {u[1] == sol[1], u'[1] == sol'[1]},
     "ImplicitSolver" -> {"Newton", 
       "StepControl" -> {"LineSearch", 

         "MaxRelativeStepSize" -> 1/110}}}];

ListLinePlot@%

I just kept decreasing "MaxRelativeStepSize" until I found something that almost worked.

Finally, sometimes one needs more control, in which case one can set up one's own shooting method with ParametricNDSolve[]. See for instance How to avoid NDSolve::ndsz problem (singularity problem)