# Solving a non-linear ODE system with ParametricNDSolveValue

I'm trying to solve a non linear ODE numerically with ParametricNDSolve, but as far as I got is shown below. My problem is to set the find root correctly. What I know is this: x' == 0, x[R] == 0, x'[R] == 0. Any help? Here is my code:

c = -0.7177;
r1 = 0.8;
r2 = 125;
R = 1.39;
f[r_] := Piecewise[{{0, 0 <= r <= r1}, {900/(1 - r1^3),
r1 < r <= 1}, {0, 1 < r <= R}}]

ps = ParametricNDSolveValue[{x''[r] + (1/r) x'[r] ==
c n Exp[-x[r]] + f[r], x' == 0, x == x0}, {x, x'}, {r, 0,
R}, {x0,n}, Method -> "StiffnessSwitching"]

ff = FindRoot[{Last[ps[x0,n]][R] == 0, First[ps[x0,n]][R] == 0}, {x0, -2}]

• How do you want to combine these three solutions of three different equations? Or is this one solution of an equation with discontinuous coefficients? Sep 28, 2018 at 7:27
• It is one equations has three regions see the image. Sep 28, 2018 at 7:33
• And why in conditions $R<b$ Sep 28, 2018 at 7:37
• it is R>b. I fixed it Sep 28, 2018 at 7:50

This problem has a solution. It is given below

c = 0.72;
h = 300;
a = 15;
b = 17;
R = 25;
f[r_] := Piecewise[{{0, 0 <= r <= a}, {(3 h)/(a^3 - b^3),
a < r <= b}, {0, b < r <= R}}]

ps = ParametricNDSolveValue[{x''[r] + 2 x'[r] ==
c n0 Exp[-x[r]] + f[r], x' == 0, x == x0},
x, {r, 0, R}, {n0, x0}];

n =
FindRoot[{ps[n0, x0][R] == 0, ps[n0, x0]'[R] == 0}, {n0, -.2}, {x0,
1}]

{n0 -> 9.19855*10^-8, x0 -> 0.585175}

{Plot[
Evaluate[Table[ps[n0, 1][r], {n0, -.2, 2, .1}]], {r, 0, R},
PlotRange -> All],
Plot[ps[n[[1, 2]], n[[2, 2]]][r], {r, 0, R}]} • Ah, now I get it. The third boundary condition is used to determine the parameter n0. Good job (and of course +1)! Sep 28, 2018 at 11:44

As Alex Trounev said, this is a second-order ODE with discontinuous right-hand side. You can use Piecewise to set up the forcing term:

rhs = Piecewise[{
{c n0 Exp[-x[r]] + (3 h)/(a^3 - b^3), a <= r < b}
},
c n0 Exp[-x[r]]
]


$$\begin{cases} \frac{3 h}{a^3-b^3}+c \,\text{n0}\, e^{-x(r)} & a\leq r

Notice that for convenience, I used c n0 Exp[-x[r]] as the default term. The full equation can be set up as

c = 0.72;
h = 300;
a = 15;
b = 17;
R = 25;
ϵ = \$MachineEpsilon;

ps = ParametricNDSolveValue[
{
x''[r] + 2 x'[r] == rhs,
x[R] == 0,
x'[R] == 0
},
x,
{r, ϵ, R},
{n0},
Method -> "StiffnessSwitching", WorkingPrecision -> 30
];


Solving it for a given parameter

f = ps[0.00001];


Plotting the result:

Plot[f[r], {r, ϵ, R}] Something must be wrong in your model: The solution blows up heavily towards $$r = 0$$ so that x' == 0 cannot be expected. Actually, you cannot prescribe more than two boundary conditions for an ODE of order 2.

• "Actually, you cannot prescribe more than two boundary conditions for an ODE of order 2." - I'm actually surprised at how frequently people don't remember to check that the number of their conditions matches the order of their ODE... Sep 28, 2018 at 8:29
• @Henrik Schumacher Look at the solution from the other side at x=0 Sep 28, 2018 at 11:41
• @user60416 What is the point of updating the code if the problem has already been solved? Do you want to solve a problem for other data? Oct 10, 2018 at 4:33
• @user60416 It is indeed not fair to invalidate already given answers by changing the question. Oct 10, 2018 at 6:40
• @user60416 You can ask a new question, and leave this one unchanged. Oct 10, 2018 at 10:31