# Solving a nonlinear systems of coupled differential equations with boundary conditions

I am trying to solve the following systems of coupled differential equations with boundary conditions (BC) at $0$ and at $∞$ :

$y_{1}'(x)=\frac{-\sqrt{\frac{2}{\pi}}\,\frac{\alpha}{18}\,x^4+x\,y_{2}(x)}{\frac{\alpha}{6}+\frac{x^3}{3}y_{4}(x)+\frac{1}{3}y_{3}(x)+\frac{2}{5}y_{5}'(x)+\frac{6}{5\,x}y_{5}(x)}y_{1}(x)$
BC : $y_{1}(0)=1 \;;\; y_{1}(\infty)=0$

$y_{2}'(x)=\sqrt{\frac{2}{\pi}}\,x^2\,y_{1}(x)$
BC : $y_{2}(0)=0\; ; \;y_{2}(\infty)=1$

$y_{3}'(x)=\sqrt{\frac{2}{\pi}}\,x^4\,y_{1}(x)$
BC : $y_{3}(0)=0\; ;\; y_{3}(\infty)=3$

$y_{4}'(x)=-\sqrt{\frac{2}{\pi}}\,x\,y_{1}(x)$
BC : $y_{4}(\infty)=0$

$\frac{\alpha}{4} \left[\frac{4}{3} x \frac{\partial}{\partial x}(\frac{1}{x^4}\frac{\partial y_{1}(x)}{\partial x})+\frac{2}{5}\frac{1}{x^7}\frac{\partial }{\partial x}\left(x^3 y_{5}(x)\right)\right]+\frac{27}{35}\frac{\alpha}{x^4}\frac{\partial}{\partial x}\left[x^3 \frac{\partial}{\partial x}(\frac{y_{5}(x)}{x^2})\right]-\frac{3}{x^3}y_{5}(x)=0$
BC : $y_{5}(0)=0\; ;\; y_{5}(\infty)=0$

Where $\alpha$ are know constants : $\alpha=0.001$

The mathematica code is:

NDSolve[{y1'[x]==((-Sqrt[(2/Pi)] alpha/18 x^4+x y2[x])/(alpha/6+6/5 1/x y5[x]+2/5 y5'[x]+ y3[x]/3+( x^3 y4[x])/3))y1[x],
y2'[x]==Sqrt[2/Pi] x^2 y1[x],
y3'[x]==Sqrt[2/Pi] x^4 y1[x],
y4'[x]==-Sqrt[(2/Pi)]x y1[x],
-(3/x^3)y5[x]+alpha/4 (4/ 3 x D[1/ (x^4) y1'[x],x]+2/ 5 1/ (x^7) D[x^3 y5[x],x])+alpha/4 36/ 7 1/ (x^4) D[3/ 5 x^3 D[1/ (x^2) y5[x],x],x]==0,y1[0]==1,y2[Infinity]==1,y3[Infinity]==3,y4[Infinity]==0,y5[Infinity]==0},{y1,y2,y3,y4,y5},{x,0,10},Method->{"Shooting","StartingInitialConditions"->{y1[0]==1,y2[Infinity]==1,y3[Infinity]==3,y4[Infinity]==0,y5[Infinity]==0}}]


I get the following output:

NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations. >>

NDSolve::ndsv: Cannot find starting value for the variable y1^\[Prime]. >>


Please, How can I find solutions of these differential equations?

Some progress can be made by noting that NDSolve cannot handle undefined alpha and boundary conditions at infinity, and that y5 appears nowhere in the first four ODEs and therefore can be decoupled from the other equations. It also makes sense to shoot from large x, because three of the four boundary conditions are there. So, we might try

xmax = 2.1;
sol = NDSolve[{y1'[x] == ((-Sqrt[(2/π)] alpha/18 x^4 + x y2[x])/
(alpha/6 + y3[x]/3 + (x^3 y4[x])/3)) y1[x],
y2'[x] == Sqrt[2/π] x^2 y1[x], y3'[x] == Sqrt[2/π] x^4 y1[x],
y4'[x] == -Sqrt[(2/π)] x y1[x], y1[0] == 1, y2[xmax] == 1,
y3[xmax] == 3, y4[xmax] == 0} /. alpha -> 10^-3,
{y1, y2, y3, y4}, {x, 0, xmax}, Method -> {"Shooting",
"StartingInitialConditions" -> {y1[xmax] ==.68477,
y2[xmax] == 1, y3[xmax] == 3, y4[xmax] == 0}}];


which produces the solution,

Plot[Evaluate[{y1[x], y2[x], y3[x], y4[x]} /. sol], {x, 0, xmax},
PlotRange -> All, AxesLabel -> {x, y}]


Unfortunately, y1[xmax] is not close to zero, suggesting that xmax has been chosen as too small a number.

So, what value of xmax is appropriate? To find out, solve for and plot the asymptotic solutions.

asym = FullSimplify[First@DSolve[{y1'[x] == ((-Sqrt[(2/π)] alpha/18 x^4)/(alpha/6)) y1[x],
y2'[x] == Sqrt[2/π] x^2 y1[x], y3'[x] == Sqrt[2/π] x^4 y1[x],
y4'[x] == -Sqrt[(2/π)] x y1[x]}, {y1[x], y2[x], y3[x], y4[x]},x]
/. {C[2] -> 1, C[3] -> 3, C[4] -> 0}, x > 0];
Plot[Evaluate[{y1[x], y2[x], y3[x], y4[x]} //.
Flatten[{C[1] -> 1, asym}]], {x, 0, 4}, AxesLabel -> {x, y}]


(C[1] -> 1 is an arbitrary choice.) To see where the asymptotic result is accurate, plot the error in the original first equation, when the asymptotic solutions are inserted.

Plot[Evaluate[(y1'[x] - ((-Sqrt[(2/π)] alpha/18 x^4 + x y2[x])/(alpha/6 + y3[x]/3 +
(x^3 y4[x])/3)) y1[x]) //. Flatten[{alpha -> 10^-3, C[1] -> 1, asym,
y1'[x] -> D[First@asym[[1,2]], x]}] // Simplify], {x, 0, 4}, AxesLabel -> {x, err}]


We see, therefore, that xmax = 3 would be a good choice for NDSolve. Unfortunately, I have not been able to find a good Shooting guess to obtain numerical solutions for xmax > 2.1. Thus, the plots above provide a qualitative picture of the solution, but not a quantitative one.

• Thank you @bbgodfrey. But I did not understand why at x=xmax, y1[xmax] == 0.68477 and not y1[xmax] ==0 as mentioned in the BC of y1. Namely, y1 is a Maxwellian solution. – Betatron Feb 27 '16 at 14:17
• @Betatron y1[xmax] is not close to zero, because xmax is too small. Unfortunately, I have not been able to obtain a solution for the more appropriate xmax = 3. The challenge is handling the possibility that alpha/6 + y3[x]/3 + x^3 y4[x]/3 may vanish somewhere. I continue to think about this issue. By the way, your question states that both y2 and y3 vanish at the origin. Are you confident that using those boundary conditions instead of the corresponding ones at large x would yield the same answer? – bbgodfrey Feb 27 '16 at 14:32
• Sorry @bbgodfrey. I made a mistake in the first differential equation. It's missing y5. Normally : y1'[x]==((-Sqrt[(2/Pi)] alpha/18 x^4+x y2[x] )/(alpha/6+6/5 1/x y5[x]+2/5 y5'[x]+ y3[x]/3+( x^3 y4[x])/3)) y1[x] where y5 satisfy the the fifth equation : -(3/ (x^3))y5[x]+alpha/4 (4/ 3 x D[1/ (x^4) y1'[x],x]+2/ 5 1/ (x^7) D[x^3 y5[x],x])+alpha/4 36/ 7 1/ (x^4) D[3/ 5 x^3 D[1/ (x^2) y5[x],x],x]==0. Thank you. Again, I apologize. – Betatron Feb 27 '16 at 14:54
• @Betatron To make any headway, you need to change dependent variables, so that only one second-derivative appears in the fifth ODE. In addition, you need to define alpha in the code and replace Infinity by 10 (or whatever upper bound of integration you choose). Even then, the first ODE creates great problems. – bbgodfrey Feb 27 '16 at 20:18
• @Betatron Have you made any progress on this question? – bbgodfrey May 18 '17 at 12:09