Solving a nonlinear systems of coupled differential equations with boundary conditions

Question

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?

Answer 1

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.