# Solving parametric ODEs at large parameter values

Question 101668 was not particularly clear and was put on hold for that reason. The author did clarify it a bit, but it seems unlikely that the question will be reopened. Through an exchange of comments with the author, I believe that I understand what he asking, which I present here and for which I offer an answer. Perhaps others can offer better answers.

The OP wishes to solve the system of equations

{f'''[x] + f[x] f''[x] + 4 - f'[x] f'[x] == 0, g''[x] + p[[i]] f[x] g'[x] == 0,
f == f' == 0, g' == -1, f' == 2, g == 0}


for a set of parameter values.

p = {0.01, 0.02, 0.04, 0.06, 0.08, 0.1, 0.2, 0.4, 0.6, 0.72, 1, 5, 7, 10, 100, 1000};


The obvious approach is

a = ParametricNDSolve[{f'''[x] + f[x] f''[x] + 4 - f'[x] f'[x] == 0,
g''[x] + p0 f[x] g'[x] == 0, f == f' == 0, g' == -1,
f' == 2, g == 0}, {f, g}, {x, 0, 5}, {p0}];


over the domain {x, 0, 5}, which works well for p as large as 7

Plot[Evaluate[{g[x] /. a, f[x] /. a}], {x, 0, 5}, PlotRange -> All] but appears to become unstable for larger values of p. (According to the OP, Maple can solve this problem for large p without difficulty.) One might hope that adding the options

Method -> "StiffnessSwitching", MaxSteps -> Infinity


would provide large p answers, but it is exceedingly slow. Setting

WorkingPrecision-> 30


instead helps a bit, but not much. It too is slow. The question, then, is - how to obtain solutions for large p in a reasonable amount of time.

One approach to obtaining large p results involves two observations. First, f is independent of p and can be solved at the outset.

sf = First[f /. NDSolve[{f'''[x] + f[x] f''[x] + 4 - f'[x] f'[x] == 0,
f == f' == 0, f' == 2}, f, {x, 0, 5}]];


Second, a first integral exists for the ODE involving g.

lgp = -Integrate[sf[x], x];


where lgp is Log[g'[x]] for p = 1. Thus, two applications of Integrate produce a solution.

hp[p0_?NumericQ] := FunctionInterpolation[-Exp[p0 lgp], {x, 0, 5},
InterpolationPoints -> 1000, InterpolationOrder -> 2]
h[p0_?NumericQ, x_?NumericQ] := -Integrate[hp[p0][x0], {x0, x, 5}]


Note that InterpolationPoints must be fairly large to provide reasonable accuracy at large p.

Plot[Evaluate[h[1000, x]], {x, 0, 5}, PlotRange -> All, AxesLabel -> {x, g}] This problem should be related to this one. Based on the experience got in that post, first analyse the equation a bit:

DSolve[{g''[x] + p0 f[x] g'[x] == 0(*, g == 0*), g' == 0}, g[x], x]

{{g[x] -> C}}


Apparently the solution will be {{g[x] -> 0}} if g == 0 (though DSolve will fail if this b.c. is directly included in it), which suggests that we can force g[x] and g'[x] to be exact 0 once they're close enough to 0 to avoid the error accumulation. So here comes the solution:

sf = NDSolveValue[{f'''[x] + f[x] f''[x] + 4 - f'[x] f'[x] == 0, f == f' == 0,
f' == 2}, f, {x, 0, 5}];

sg = With[{p0 = 1000, f = sf, e = \$MachineEpsilon},
NDSolveValue[{g''[x] + p0 f[x] g'[x] == 0, g' == -1, g == 0,
WhenEvent[{g[x] < e, g'[x] < e}, {g[x] -> 0, g'[x] -> 0}]}, g, {x, 0, 5}]]

Plot[sg[x], {x, 0, 5}] 