# Solve a Second-Order Differential Equation Numerically, with boundary conditions?

Below is an ODE with BC define as x[R]=0, x'[0]=0, x'[R]=0 and parameter n. The ODE is stiff at certain data, and I need to see the behavior of x' and x for given parameters. I appreciate any help.

c = -1;
r1 = 0.8;
r2 = 1;
R = 1.29;
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] == 0, x[0] == 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}]

• I have not yet found a convergent solution. It is possible that it is necessary to change the boundary conditions. – Alex Trounev Oct 10 '18 at 18:10
• @AlexTrounev see my update – user60416 Oct 10 '18 at 19:09
• The figure shows R>r2, and in the data otherwise r2=125>R=1.29. – Alex Trounev Oct 11 '18 at 3:53
• Here r1/r2=0.8 and R/r2 =1.29, so instead the integral from r1<r<r2 is it from r1<r<1 where 1 = r2. – user60416 Oct 11 '18 at 4:40
• What is the meaning of the parameter с, if the solution of the equation depends on the product c * n? – Alex Trounev Oct 11 '18 at 4:47

The solution of the equation depends on the product c * n, but not on c. Therefore, we can set c = -1, since the parameter n is found in the solution process.

c = -1;
r1 = 8/10;
r2 = 1;
R = 129/100;
r0 = 10^-6;
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'[r0] == 0, x[r0] == x0},
x, {r, r0, R}, {x0, n}, Method -> "StiffnessSwitching",
WorkingPrecision -> 30];
ps1 = ParametricNDSolveValue[{x''[r] + (1/r) x'[r] ==
c n Exp[-x[r]] + f[r], x'[r0] == 0, x[r0] == x0},
x', {r, r0, R}, {x0, n}, Method -> "StiffnessSwitching",
WorkingPrecision -> 30]

ff = FindRoot[{ps[x0, n][R] == 0, ps1[x0, n][R] == 0}, {x0,
1.27717130}, {n, 393/10}]

{x0 -> 1.27686, n -> 39.2879}

{Plot[Evaluate[ps[Last[ff[[1]]], Last[ff[[2]]]][r]], {r, r0, R},
PlotRange -> All, AxesLabel -> {"r", "x"}],
Plot[Evaluate[ps1[Last[ff[[1]]], Last[ff[[2]]]][r]], {r, r0, R},
PlotRange -> All, AxesLabel -> {"r", "x'"}]}


• How do you estimate theses numbers {x0, 1.27717130}, {n, 393/10}]? – user60416 Oct 11 '18 at 19:23
• Would you please explain to me this line Evaluate[ps[Last[ff[[1]]], Last[ff[[2]]]][r]]? like what is 2 refer to? – user60416 Oct 12 '18 at 1:13
• To find the initial values {x0,n}, I used an iteration prediction procedure. The second question is not clear to me. Why don't you just see what corresponds to what? ps->x,ps1->x',Last[ff[[1]]]-> 1.27686,Last[ff[[2]]]->39.2879 – Alex Trounev Oct 12 '18 at 2:57
• Thanks! Could you share your code of (iteration prediction procedure)? – user60416 Oct 12 '18 at 3:31
• Unfortunately I can not give the code. But it is not difficult to write in the Wolfram Language. – Alex Trounev Oct 12 '18 at 3:35