Shooting Method with extra Unknown Condition

I've been trying to solve a system of 3 coupled 2nd order ODEs, for a real variable $$x$$, $$0\geq x\leq \infty$$. The equations are the following: \begin{align} &x^{2}\,h''(x) - x\,h'(x) + x^{2}\,g^{2}(x)\left[1-h(x)\right] = 0\,,\\ &x^{2}\,f''(x) + x\,f'(x) - \lambda\, x^{2}\,f(x)\left[f^{2}(x) + g^{2}(x) - 2\right] = 0\,,\\ &x^{2}\,g''(x) + x\,g'(x) - \frac{1}{2}\,g(x)\left[1-h(x)\right]^{2} - \lambda\, x^{2}g(x)\left[f^{2}(x) + g^{2}(x) - 2\right] = 0\,. \end{align}

In addition, the BCs are (where my problem starts): $$h(0)=0=g(0)\,, \quad f(0)=\Omega$$ and $$h(x\to\infty)=f(x\to\infty)=g(x\to\infty)=1\,.$$

First of all, I decided to solve to some finite $$x$$ such as $$x_{max}$$ and then try to increase this domain. Then, my problem consists of how I should "tell" NDSolve that I don't know the value $$\Omega$$ is going to have. I know that, somehow, the numerical solution must find an appropriate value for $$\Omega$$ that agrees with the whole solution. But I cannot understand how I can do this. My starting code is the following:

lambda = 0.5;
eps = 0.001;
xmax = 5;
eq1=x^2*h''[x] - x*h'[x] + x^2*(g[x]^2) (1 - h[x]);
eq2= x^2*f''[x] + x*f'[x] - lambda*x^2*f[x] ((f[x]^2) + (g[x]^2) - 2);
eq3= x^2*g''[x] + x*g'[x] -  1/2*g[x] (1 - h[x])^2 - lambda*x^2*g[x] ((f[x]^2) + (g[x]^2) -
2);

sols=First[NDSolve[{eq1==0,eq2==0, eq3== 0,h[eps] == 0,f[eps] == Omega, g[eps]==0}, {f[x],
g[x], h[x]}, {x, eps, xmax},Method -> {"Shooting","StartingInitialConditions" -> {h[eps] ==
0,f[eps] == Omega, g[eps] == 0}}, WorkingPrecision -> 5]];


As you can see, my code is incomplete. The shooting method would need 6 initial conditions for the (converted) IVP -> 3 from the BC at $$x=0$$ and the shooting for the 3 first-order derivatives. However, since I don't know (a priori) the value of $$\Omega$$, I'm stuck :(

Ps.: From my problem, I know I can put by hand that all the first-order derivatives go to zero when $$x\to\infty$$... But this would add too many conditions for Mathematica, right?

Could you, please, provide any advice on how I can tackle the problem?

• For problems like these, it might be better to use a mesh and apply Newton's cooling method starting from an initial guess for the solution, rather than using the shooting method. Commented Sep 18, 2020 at 22:24
• Hi @QuantumDot. I do not have experience with numerical solutions. Would this "Newton's cooling method" similar to what Numerical Recipes call Relaxation? And I have tried to solve this by shooting because the paper I was reading told the readers those equations were solved by a shooting method together with a 5th order RK for integration. However, I can try to learn relaxation, too. Thanks :) Commented Sep 19, 2020 at 23:29

With NDSolve and Method -> "Shooting", I was unable to integrate past x = 2.93. With Method -> "FiniteElement", I was able to integrate as far as x = 7.5, but omega could only be approximated. In general, the problem appeared to be inadequate numerical accuracy near x = 0. The following worked much better.

First, obtain approximate symbolic solutions near x = 0

sh0 = DSolveValue[{eq1 == 0, h[0] == 0} /. {f[x]^2 -> omega^2, g[x]^2 -> 0}, h[x], x] /.
C[1] -> ch
(* (ch x^2)/2 *)

sf0 = Simplify[DSolveValue[{eq2 == 0, f[0] == omega} /. {f[x]^2 -> omega^2, g[x]^2 -> 0},
f[x], x], omega^2 < 2]
(* omega BesselJ[0, (Sqrt[2 - omega^2] x)/Sqrt[2]] *)

sg0 = Simplify[DSolveValue[{eq3 == 0, g[0] == 0} /. {f[x]^2 -> omega^2,
g[x]^2 -> 0, h[x] -> 0}, g[x], x], omega^2 < 2] /. C[1] -> cg
(* cg (BesselJ[1/Sqrt[2], (Sqrt[2 - omega^2] x)/Sqrt[2]] - (BesselJ[1/Sqrt[2], 0]
BesselY[1/Sqrt[2], (Sqrt[2 - omega^2] x)/Sqrt[2]])/BesselY[1/Sqrt[2], 0]) *)


Then, obtain the three parameters {omega, cg, ch} by what might be called do-it-yourself shooting, with the symbolic solutions near x = 0 as initial conditions.

xmax = 10; eps = .1;
sp = ParametricNDSolveValue[{eq1 == 0, eq2 == 0, eq3 == 0,
{f[x] == sf0, g[x] == sg0, h[x] == sh0, f'[x] == D[sf0, x], g'[x] == D[sg0, x],
h'[x] == D[sh0, x]} /. x -> eps}, {f[xmax], g[xmax], h[xmax]}, {x, eps, xmax},
{omega, cg, ch}, PrecisionGoal -> 10, AccuracyGoal -> 10];

FindRoot[sp[omega, cg, ch] - 1, {{omega, 1.35617}, {cg, 1.3415}, {ch, 0.325786}},
Evaluated -> False]
(* {omega -> 1.35617, cg -> 1.3415, ch -> 0.325786} *)


Finally, compute and solve the equations with the parameters just determined.

NDSolveValue[{eq1 == 0, eq2 == 0, eq3 == 0, {f[x] == sf0, g[x] == sg0, h[x] == sh0,
f'[x] == D[sf0, x], g'[x] == D[sg0, x], h'[x] == D[sh0, x]} /. x -> eps} /. %,
{f[x], g[x], h[x]}, {x, eps, xmax}];
Plot[%, {x, eps, xmax}, ImageSize -> Large, AxesLabel -> {x, "f,g,h"},
LabelStyle -> {15, Bold, Black}]


The initial guesses for FindRoot were obtained by integrating the equations for xmax = 3 and using the result as guesses for xmax = 4, etc. Note that even xmax = 10 is not in the asymptotic domain of the equations. Increasing xmax to, say 20 undoubtedly would require higher WorkingPrecision and an automated process for gradually increasing xmax, both of which are feasible.

• Thank you so much! It has helped me a lot. Commented Sep 24, 2020 at 14:01
• @bbgodfrey The term sp[omega, cg, ch] - 1 in FindRoot, why is there a -1? Also, supposedly you didn't have much idea on the problem how did you choose your starting point for FindRoot? I have tried changing the initial point a bit and it produced errors like ParametricNDSolveValue::ndsz InterpolatingFunction::dmval FindRoot::lstol Commented Nov 25, 2022 at 16:52
• @mathemania 1 is the desired asymptotic value for the three dependent variables. I used FindRoot to obtain initial conditions that achieved this value for my solution. The final paragraph of my answer explains how I obtained my initial guesses. Commented Nov 25, 2022 at 18:10