# Nonlinear system of ODEs with boundary conditions

I'm trying solve this problem:

g'(r) = a(r)g(r)/r, (1/r)a'(r) = g(r)^2-1

which have the following boundary conditions: a(0)=n and g(infinity)=1, where n>=1.

I have problem with the error for n=1[hence, failed for the other n values as well.]:

"Infinite expression 1/0.^2 encountered"

using DSolve directly.

A work colleague solved for the values n=1,2,5,10 using the software MAPLE. But, I don't know to use MAPLE and I would like to get these results by Mathematica. Below are the plots found by MAPLE. • Please include the Mathematica code text in your question. – xzczd Aug 8 '18 at 16:58
• Welcome to Mma.SE. Start by taking the tour now and learning about asking and what's on-topic. Always edit if improvable, show due diligence, give brief context, include minimal working example of code and data you are working on in formatted form. By doing all this you help us to help you and likely you will inspire great answers. The site depends on participation, as you receive give back: vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned – rhermans Aug 8 '18 at 17:10
• See this previous answer. Compactifying in MMA is rather trivial. – Hector Aug 8 '18 at 17:10
• I find that DSolve does not produce an error message but instead returns unevaluated, meaning that it cannot solve the problem.. – bbgodfrey Aug 8 '18 at 21:34
• Maple code:sol := proc (n) options operator, arrow; dsolve([(D(g))(r) = a(r)*g(r)/r, (D(a))(r)/r = g(r)^2-1, a(1/100000) = n, g(10) = 1], numeric) end proc;plots:-display(seq({plots:-odeplot(sol(n), [r, g(r)], linestyle = n, color = ColorTools:-Color([1/n, 1/(1+n), 1/(n+4)]))}, n = 1 .. 7));plots:-display(seq({plots:-odeplot(sol(n), [r, a(r)], linestyle = n, color = ColorTools:-Color([1/n, 1/(1+n), 1/(n+4)]))}, n = 1 .. 7)); – Mariusz Iwaniuk Aug 11 '18 at 17:17

Another way:

eps = 10^-10; end = 10;
sol[n_] := NDSolve[
{g'[r] == a[r]*g[r]/r, a'[r]/r == g[r]^2 - 1, a[eps] == n, g[end] == 1}
,{a[r], g[r]}, {r, eps, end}, Method -> {"Shooting",
"StartingInitialConditions" -> {a[end] == eps, g[end] == (1 - eps)}}];

Plot[Evaluate[g[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end},PlotRange -> All,
PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, PlotLabels -> "g(r)"]
Plot[Evaluate[a[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end}, PlotRange -> All,
PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, PlotLabels -> "a(r)"] For larger values e.g.end=20:

eps = 10^-20; end = 20;
sol[n_] := NDSolve["Code the same",WorkingPrecision -> 25];
Plot["Code the same"] // Quiet(*Error messages form InterpolatingFunction to Quiet *)

• Thanks for your answer, you've been a big help. – sotnasac Aug 12 '18 at 18:49
• Dear @MariuszIwaniuk, I was very grateful for the help on this occasion. Now, if it is possible to help again, how could I make a projection of these solutions in the plan? That is, would it be possible to make a density plot? – sotnasac Jul 19 '19 at 19:19

Error in boundary condition corrected

To solve this ODE system, begin by determining its behavior near r == 0.

Series[{g'[r] - (a[r] g[r])/r, a'[r]/r + 1 - g[r]^2}, {r, 0, 2}] // Normal;
Thread[Flatten@CoefficientList[r % /. a -> n, r] == 0];
Solve[%, {a', a'', a''', a'''', g, g', g'', g'''}]

(* a' -> 0, a'' -> -1, a''' -> 0, a'''' -> 0,
g -> 0, g' -> 0, g'' -> 0, g''' -> 0  *)


In other words, g[r] vanishes to all orders at the origin for most values of n. However, if n is a positive integer, the same set of commands leaves Derivative[n][g] undetermined. So, g[eps] is approximately proportional to eps^n. Next, solve the equations numerically by shooting, using the method previously employed here and elsewhere on this site. (As noted in my comment above, DSolve returns unevaluated when applied to these equations.) A typical calculation for n == 5 is given by

eps = 10^-5; end = 20;
sp = ParametricNDSolveValue[{g'[r] == a[r] g[r]/r, a'[r]/r == g[r]^2 - 1,
a[eps] == n0, g[eps] == eps^n0 gp,
WhenEvent[g[r] > 101/100, {bool = 1, "StopIntegration"}],
WhenEvent[a[r] < 0, {bool = 0, "StopIntegration"}]},
{a[r], g[r]}, {r, eps, end + 5}, {gp, wp0, n0},
WorkingPrecision -> wp0, Method -> "StiffnessSwitching"];
bl = 0; bu = 2; imax = 100; wp = 75; n = 5;
Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp, n];
rm = (First[st] /. r -> "Domain")[[1, 2]];
If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i;
If[bool == -1, Return[]], {i, imax}];
s[n] = st;


For larger values of n, wp must be increased. For instance, with n = 10, wp = 90 is required. Results for n = {1, 2, 5, 10} are plotted below.

Plot[{s, s, s, s}, {r, eps, 10}, PlotRange -> All, AxesLabel -> {r, "a, g"},
ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}] The solutions takes no more than several seconds each.

• I have appreciated your answer. Thank you very much @bbgodfrey. You've already been very helpful. A friend got solve also for values n>1, but he used the software MAPLE. I don't know to use the MAPLE, but I'll try to use your answer like a procedure. – sotnasac Aug 10 '18 at 1:49
• @AndréCavalcante Because the ODEs are nonlinear, not all boundary conditions are consistent with them. In this case, I believe that A == 1 is the only consistent boundary condition, unless you are willing to accept g[r] == 0 as a solution. Perhaps, I could get Mathematica to produce numerical solutions for other values of n, but they would be incorrect, I believe. – bbgodfrey Aug 10 '18 at 1:58
• I have edited my question and added the plots obtained by MAPLE. Indeed, using its answer it was not possible to solve for values greater than n == 1. In addition, by behavior near r == 0, I have noticed that for n == 2, for example, leaves g'' undetermined. – sotnasac Aug 11 '18 at 16:58
• @AndréCavalcante I just realized that I had made an error in the g[eps] boundary condition. In fact, any positive integer value of n gives a valid solution. Sorry. – bbgodfrey Aug 11 '18 at 18:17
• Again, thanks for your help. – sotnasac Aug 12 '18 at 18:42