# Need help solving a system of two 1st order nonlinear differential equations

The original system of equations reads:

$\begin{cases} f'(r) + f(r) \left(a(r) - \frac{1}{r}\right) = 0,\\ f^2(r) + a'(r) + \frac{a(r)}{r} - 1 = 0\,, \end{cases}$

with boundary conditions $f(0) = 0\,,\,\, f(\infty) = 1\,.$ I know that it allows non-divergent nontrivial solutions $\left(0 \leq f(r) \leq 1\,,\,\, f(r \to 0) \sim r\right.,$ and approaching exponentially to $1$ at large $r$$\left. \right)$.

Due to the boundary condition at infinity I tried the following change of variables: $r = \frac{1}{x} - 1$, which transforms it to the following system:

$\begin{cases} -x^2 f'(x) + f(x) \left(a(x) - \frac{x}{1 - x}\right) = 0,\\ f^2(x) - x^2 a'(x) + a(x) \frac{x}{1 - x} - 1 = 0\,, \end{cases}$

with boundary conditions $f(0) = 1\,,\,\, f(1) = 0\,,$ which Matematica (NDSolve) can't solve:

NDSolve[{-x^2 f'[x] + f[x] (a[x] - x/(1 - x)) == 0, f[x]^2 - x^2 a'[x] + a[x] x/(1 - x) == 1,
f[0] == 1, f[1] == 0}, {f,a}, {x, 0, 1}]


One can eliminate $a(x)$ to obtain a 2nd-order differential equation for $f(x)$:

$f^2(x) + x^4\left( \frac{f'(x)}{f(x)}\right)^2 - x^4 \frac{f''(x)}{f(x)} +x^3\left( \frac{1}{1 - x} - 2\right)\frac{f'(x)}{f(x)} - 1 = 0\,,$ for which NDSolve again has no answer:

NDSolve[{f[x]^2 + x^4 (f'[x]/f[x])^2 - x^4 f''[x]/f[x] + x^3 (1/(1 - x) - 2) f'[x]/f[x] - 1 == 0,
f[0] == 1, f[1] == 0}, f, {r, 0, 1}]


Any suggestions how to solve this problem in Mathematica (and not Matlab)?

There is a problem with boundary conditions. Changing them and fixing a few typos you can get:

eq = {-f[x]^2 + f[x]^4 +
x^4 Derivative[1][f][x]^2 + (x^3) /(-1 + x)
f[x] ((1 - 2 x) f'[x] - (-1 + x) x f''[x]) == 0, f[1/10] == 1,
f'[1/10] == -1/10};


s = NDSolve[eq, f, {x, 1/10, 9/10}];
Plot[Evaluate[f[x] /. s], {x, 1/10, 9/10}, PlotRange -> All]


Usually people first try to consider some expansions, approximations, etc. at the boundaries to understand the behaviour. So the system solvable - just understand it first and set up proper BCs. Perhaps there is an error somewhere in derivation of final equation that it stops corresponding to your BCs?

• I don't believe there's a mistake. You are right, however, that boundary conditions are very important here. The solution I was after looks different from yours because its shape is highly dependent on BCs. Note that as you decrease the value of the first derivative while the boundary slowly approaches 0 from 1/10, the shape of the solution changes drastically. It does not even seem to converge when doing it this way. Instead, you can get the solution found in literature by assuming the following BCs: f[1 - 10^-5] == 0.85329*10^-5, f'[1 - 10^-5] == -0.85329. Why is it so unstable/sensitive? – Alexey Galda Aug 19 '14 at 3:39

Starting from the two equations in the question,

eq1 = D[f[r], r] + f[r] (a[r] - 1/r);
eq2 = f[r]^2 + D[a[r], r] + a[r]/r - 1;


eliminate a[r].

arule = Solve[eq1 == 0, a[r]][[1, 1]];
darule = D[%, r] // Simplify;
eq3 = -f[r] (eq2 /. {arule, darule}) // Simplify
(* f[r] - f[r]^3 + Derivative[1][f][r]/r - Derivative[1][f][r]^2/f[r] +
Derivative[2][f][r] *)


Next, formulate the r = 0 boundary condition by expanding f[r] as a power series there and requiring that f[0] == 0.

Series[eq3 /. Thread[Rule[#, # /. f -> Function[{r}, Sum[c[i] r^i, {i, 4}]]] & /@
{f[r], Derivative[1][f][r], Derivative[2][f][r]}], {r, 0, 2}] // Normal
(* c[2] + r (c[1] - c[2]^2/c[1] + 4 c[3]) +
r^2 (c[2] + c[2]^3/c[1]^2 - (4 c[2] c[3])/c[1] + 9 c[4]) *)


Evidently, c[1] is arbitrary, while c[2] (and c[4]) is zero.

From question 91854 there are a variety of related ways to proceed. A particularly simple approach, due to Michael E2, is as follows.

r0 = 10^-5; r1 = 50;
sEvent = ParametricNDSolveValue[{eq3 == 0, f[r0] == f'[r0] r0,
f'[r0] == b, WhenEvent[f'[r] < 0, "StopIntegration"],
WhenEvent[f[r] > 1, "StopIntegration"]}, f, {r, r0, r1}, {b}, WorkingPrecision -> 45];
dist[b_?NumericQ] := sEvent[b]["Domain"][[1, 2]]
FindMaximum[dist[b], {b, 8529/10000}, WorkingPrecision -> 45];
{dmax, b0} = {First@%, b /. Last[%]}
(* {19.7939826947595284038626883012992096669595719,
0.853177865433588923966572213040425596262318851} *)

Plot[sEvent[b0][r], {r, r0, dmax}, PlotRange -> All, AxesLabel -> {r, f}]


For completeness, a[r] can be obtained from

arule[[2]] // Apart
(* 1/r - Derivative[1][f][r]/f[r] *)

Plot[1/r - sEvent[b0]'[r]/sEvent[b0][r], {r, r0, dmax}, AxesLabel -> {r, a}]