# Solution of Coupled second-order ODEs and plot the diagram

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 == 10^-8},
{y2'[x] == y4[x], y2 == 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 == 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}
}, {y1[x], y2[x], y3[x], y4[x]}, x]


• This doesn't look like something DSolve[] can deal with; why not use NDSolve[] instead if all you need is a plot? – J. M.'s technical difficulties May 25 '16 at 5:52
• Gahir: @J.M. wasn't asking for clarification. He was suggesting that you try to use NDSolve instead of DSolve, since your equations are unlikely to be analytically solvable by Mathematica. Look up the documentation for NDSolve for the correct syntax and how to plot the outputs. – march May 25 '16 at 6:05

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 == 10^-8;
C2 = y2 == 0;
C3 = y1' == 0;
C4 = y2' == 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]