# 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. – QuantumDot Sep 18 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 :) – MLPhysics Sep 19 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. – MLPhysics Sep 24 at 14:01