Solution of Coupled second-order ODEs and plot the diagram

Question

We have two second-order Coupled differential equations as the followings:

$$\left\{\begin{array}{lr} \displaystyle \frac{{{d^2}{y_1}}}{{d{x^2}}} = \{ \frac{{\sqrt {\frac{{1 - {\varepsilon ^2}}}{{{{(\sqrt {1 + {a_0}^2} + {y_2})}^2} - {\varepsilon ^2}(1 + {y_1}^2)}}} }}{{{\varepsilon ^2}}} - {\omega ^2}\} {y_1}, \\ \displaystyle \frac{{{d^2}{y_2}}}{{d{x^2}}} = \frac{{\sqrt {\frac{{1 - {\varepsilon ^2}}}{{{{(\sqrt {1 + {a_0}^2} + {y_2})}^2} - {\varepsilon ^2}(1 + {y_1}^2)}}} \left( {\sqrt {1 + {a_0}^2} + {y_2}} \right) - 1}}{{{\varepsilon ^2}}} \end{array}\right.$$

We need to obtain $y_1$ and $y_2$ in terms of $x$ and then to draw their diagrams in terms of $x$. It should be noted that the $ω$, $α_0$ and $\varepsilon$ are some constants having the values of $1.83465945$, $0$ and $0.5$ respectively. Now, to the ease the task, we’re going to write the two aforementioned second ordered equations as four first ordered equations as the followings: \begin{align} \frac{{d{y_1}}}{{dx}} &= {y_3}, \\ \frac{{d{y_2}}}{{dx}} &= {y_4}, \\ \frac{{d{y_3}}}{{dx}} &= \{ \frac{{\sqrt {\frac{{1 - {\varepsilon ^2}}}{{{{(\sqrt {1 + {a_0}^2} + {y_2})}^2} - {\varepsilon ^2}(1 + {y_1}^2)}}} }}{{{\varepsilon ^2}}} - {\omega ^2}\} {y_1},\\ \frac{{d{y_4}}}{{dx}} & = \frac{{\sqrt {\frac{{1 - {\varepsilon ^2}}}{{{{(\sqrt {1 + {a_0}^2} + {y_2})}^2} - {\varepsilon ^2}(1 + {y_1}^2)}}} \left( {\sqrt {1 + {a_0}^2} + {y_2}} \right) - 1}}{{{\varepsilon ^2}}} \end{align} So, in this case, we need to obtain the $y_1$,$y_2$, $y_3$ and $y_4$ in terms of $x$, and only to draw $y_1$ and $y_2$ in terms of $x$. It also should be noted that the initial values of $y_1$, $y_2$, $y_3$ and $y_4$ are the followings: \begin{align} y_1(0) &= 10^{-8},\\ y_2(0) &= 0,\\ y_3(0) &= 9.81263*10^{- 10},\\ y_4(0) &= 0. \end{align}

The try I've made in Mathematica is like the following, I'm in trouble with the solution and the depiction of this problem:

epsil = 0.5;
asefr = 0;
omeg = 1.834659451;
DSolve[{
{y1'[x] == y3[x], y1[0] == 10^-8},
{y2'[x] == y4[x], y2[0] == 0},
{y3'[x] == (-omeg^2 + Sqrt[1 - epsil^2]/(epsil^2 Sqrt[(Sqrt[1 + asefr^2] + y2[x])^2 - epsil^2 (1 + (y1^2)[x])]))*y1[x], y3[0] == 9.81263*10^-10},
{y4'[x] == ((Sqrt[1 + asefr^2] + y2[x]) Sqrt[1 - epsil^2])/(epsil^2 Sqrt[(Sqrt[1 + asefr^2] + y2[x])^2 - epsil^2 (1 + (y1^2)[x])]) + 1/epsil^2, y4[0] == 0}
}, {y1[x], y2[x], y3[x], y4[x]}, x]

I do appreciate your help, in advance.

Answer 1

As mentioned in the comments that DSolve will be unable to solve this nonlinear coupled system of two ODE's. Instead you can use NDSolve. Here is my try

omega = 1.83465945;
a0 = 0;
epsilon = 0.5;

Eq1 = y1''[x] == (Sqrt[(1 - epsilon^2)/((Sqrt[1 + a0^2] + y2[x])^2 - 
           epsilon^2*(1 + y1[x]^2))]/epsilon^2 - omega^2) y1[x];
Eq2 = y2''[x] == (Sqrt[(1 - epsilon^2)/((Sqrt[1 + a0^2] + y2[x])^2 - 
           epsilon^2*(1 + y1[x]^2))]*(Sqrt[1 + a0^2] + y2[x]) - 1)/
    epsilon^2;

I choose the other two conditions randomly

C1 = y1[0] == 10^-8;
C2 = y2[0] == 0;
C3 = y1'[0] == 0;
C4 = y2'[0] == 0;

sol = NDSolve[{Eq1, Eq2, C1, C2, C3, C4}, {y1, y2}, {x, 0, 5}]


Plot[Evaluate[{y1[x] /. sol, y2[x] /. sol}], {x, 0, 5}, PlotRange -> All]