How can I solve the following higher order ODE?

Question

I am trying to solve a system of two second-order ODEs. After separating them, I obtained a fourth-order independent ODE as illustrated below. I wonder if there is a specific technique to solve it.

$$y^{(4)}+\frac{a_1}{x} y^{(3)}+\frac{a_2}{x^2}y^{(2)}+a_3y^{(2)}+\frac{a_4}{x}y^{(1)}+a_5y=0$$

Answer 1

Hint:

$y^{(4)}+\dfrac{a_1}{x}y^{(3)}+\dfrac{a_2}{x^2}y^{(2)}+a_3y^{(2)}+\dfrac{a_4}{x}y^{(1)}+a_5y=0$

$x^2y''''+a_1xy'''+(a_2+a_3x^2)y''+a_4xy'+a_5x^2y=0$

Another notable special cases appear when $a_2=0$ , the ODE reduces to

$xy''''+a_1y'''+a_3xy''+a_4y'+a_5xy=0$

which belongs to a linear ODE with linear coefficients. The solution can be found by assuming a suitable integral kernel which is so-called the “integral kernel method”.

Let $y=\int_m^ne^{xs}K(s)~ds$ ,

Then $x\int_m^ns^4e^{xs}K(s)~ds+a_1\int_m^ns^3e^{xs}K(s)~ds+a_3x\int_m^ns^2e^{xs}K(s)~ds+a_4\int_m^nse^{xs}K(s)~ds+a_5x\int_m^ne^{xs}K(s)~ds=0$

$\int_m^n(s^4+a_3s^2+a_5)e^{xs}K(s)~d(xs)+\int_m^n(a_1s^3+a_4s)e^{xs}K(s)~ds=0$

$\int_m^n(s^4+a_3s^2+a_5)K(s)~d(e^{xs})+\int_m^n(a_1s^3+a_4s)e^{xs}K(s)~ds=0$

$[(s^4+a_3s^2+a_5)e^{xs}K(s)]_m^n-\int_m^ne^{xs}~d((s^4+a_3s^2+a_5)K(s))+\int_m^n(a_1s^3+a_4s)e^{xs}K(s)~ds=0$

$[(s^4+a_3s^2+a_5)e^{xs}K(s)]_m^n-\int_m^n((s^4+a_3s^2+a_5)K'(s)+(4s^3+2a_3s)K(s))e^{xs}~ds+\int_m^n(a_1s^3+a_4s)e^{xs}K(s)~ds=0$

$[(s^4+a_3s^2+a_5)e^{xs}K(s)]_m^n-\int_m^n((s^4+a_3s^2+a_5)K'(s)+((4-a_1)s^3+(2a_3-a_4)s)K(s))e^{xs}~ds=0$

$\therefore(s^4+a_3s^2+a_5)K'(s)+((4-a_1)s^3+(2a_3-a_4)s)K(s)=0$

$(s^4+a_3s^2+a_5)K'(s)=((a_1-4)s^3+(a_4-2a_3)s)K(s)$

$\dfrac{K'(s)}{K(s)}=\dfrac{(a_1-4)s^3+(a_4-2a_3)s}{s^4+a_3s^2+a_5}$

$\int\dfrac{K'(s)}{K(s)}~ds=\int\dfrac{(a_1-4)s^3+(a_4-2a_3)s}{s^4+a_3s^2+a_5}~ds$

$\int\dfrac{d(K(s))}{K(s)}=\dfrac{1}{2}\int\dfrac{(a_1-4)s^2+a_4-2a_3}{\left(s^2+\dfrac{a_3}{2}\right)^2+a_5-\dfrac{a_3^2}{4}}~d\left(s^2+\dfrac{a_3}{2}\right)$

Case $1a$: $a_5=\dfrac{a_3^2}{4}$

Then $\int\dfrac{d(K(s))}{K(s)}=\dfrac{1}{2}\int\dfrac{(a_1-4)s^2+a_4-2a_3}{\left(s^2+\dfrac{a_3}{2}\right)^2}~d\left(s^2+\dfrac{a_3}{2}\right)$

$\int\dfrac{d(K(s))}{K(s)}=\dfrac{1}{2}\int\dfrac{(a_1-4)\left(s^2+\dfrac{a_3}{2}\right)+a_4-2a_3-\dfrac{a_3}{2}(a_1-4)}{\left(s^2+\dfrac{a_3}{2}\right)^2}~d\left(s^2+\dfrac{a_3}{2}\right)$

$\int\dfrac{d(K(s))}{K(s)}=\dfrac{1}{2}\int\left(\dfrac{a_1-4}{s^2+\dfrac{a_3}{2}}+\dfrac{a_4-\dfrac{a_1a_3}{2}}{\left(s^2+\dfrac{a_3}{2}\right)^2}\right)d\left(s^2+\dfrac{a_3}{2}\right)$

$\ln K(s)=\dfrac{a_1-4}{2}\ln\left(s^2+\dfrac{a_3}{2}\right)+\dfrac{a_1a_3-2a_4}{4s^2+2a_3}+c_1$

$K(s)=c\left(s^2+\dfrac{a_3}{2}\right)^\frac{a_1-4}{2}e^\frac{a_1a_3-2a_4}{4s^2+2a_3}$

$\therefore y=\int_m^nc\left(s^2+\dfrac{a_3}{2}\right)^\frac{a_1-4}{2}e^{xs+\frac{a_1a_3-2a_4}{4s^2+2a_3}}~ds$

But since the above procedure in fact suitable for any complex number $s$ ,

$\therefore y_p=\int_{m_p}^{n_p}c_p\left((k_pt)^2+\dfrac{a_3}{2}\right)^\frac{a_1-4}{2}e^{x(k_pt)+\frac{a_1a_3-2a_4}{4(k_pt)^2+2a_3}}~d(k_pt)=k_pc_p\int_{m_p}^{n_p}\left(k_p^2t^2+\dfrac{a_3}{2}\right)^\frac{a_1-4}{2}e^{k_pxt+\frac{a_1a_3-2a_4}{4k_p^2t^2+2a_3}}~dt$

For some $x$-independent real number choices of $m_p$ and $n_p$ and $x$-independent complex number choices of $k_p$ such that:

$\lim\limits_{t\to m_p}\left(k_p^2t^2+\dfrac{a_3}{2}\right)^\frac{a_1}{2}e^{k_pxt+\frac{a_1a_3-2a_4}{4k_p^2t^2+2a_3}}=\lim\limits_{t\to n_p}\left(k_p^2t^2+\dfrac{a_3}{2}\right)^\frac{a_1}{2}e^{k_pxt+\frac{a_1a_3-2a_4}{4k_p^2t^2+2a_3}}$

$\int_{m_p}^{n_p}\left(k_p^2t^2+\dfrac{a_3}{2}\right)^\frac{a_1-4}{2}e^{k_pxt+\frac{a_1a_3-2a_4}{4k_p^2t^2+2a_3}}~dt$ converges

Case $1b$: $a_5>\dfrac{a_3^2}{4}$

Then $\ln K(s)=\dfrac{a_1-4}{4}\ln(s^4+a_3s^2+a_5)+\dfrac{a_4-2a_3}{\sqrt{4a_5-a_3^2}}\tan^{-1}\dfrac{2s^2+a_3}{\sqrt{4a_5-a_3^2}}+c_1$

$K(s)=c(s^4+a_3s^2+a_5)^\frac{a_1-4}{4}e^{\frac{a_4-2a_3}{\sqrt{4a_5-a_3^2}}\tan^{-1}\frac{2s^2+a_3}{\sqrt{4a_5-a_3^2}}}$

$\therefore y=\int_m^nc(s^4+a_3s^2+a_5)^\frac{a_1-4}{4}e^{xs+\frac{a_4-2a_3}{\sqrt{4a_5-a_3^2}}\tan^{-1}\frac{2s^2+a_3}{\sqrt{4a_5-a_3^2}}}~ds$

But since the above procedure in fact suitable for any complex number $s$ ,

$\therefore y_p=\int_{m_p}^{n_p}c_p((k_pt)^4+a_3(k_pt)^2+a_5)^\frac{a_1-4}{4}e^{xk_pt+\frac{a_4-2a_3}{\sqrt{4a_5-a_3^2}}\tan^{-1}\frac{2(k_pt)^2+a_3}{\sqrt{4a_5-a_3^2}}}~d(k_pt)=k_pc_p\int_{m_p}^{n_p}(k_p^4t^4+a_3k_p^2t^2+a_5)^\frac{a_1-4}{4}e^{k_pxt+\frac{a_4-2a_3}{\sqrt{4a_5-a_3^2}}\tan^{-1}\frac{2k_p^2t^2+a_3}{\sqrt{4a_5-a_3^2}}}~dt$

For some $x$-independent real number choices of $m_p$ and $n_p$ and $x$-independent complex number choices of $k_p$ such that:

$\lim\limits_{t\to m_p}(k_p^4t^4+a_3k_p^2t^2+a_5)^\frac{a_1}{4}e^{k_pxt+\frac{a_4-2a_3}{\sqrt{4a_5-a_3^2}}\tan^{-1}\frac{2k_p^2t^2+a_3}{\sqrt{4a_5-a_3^2}}}=\lim\limits_{t\to n_p}(k_p^4t^4+a_3k_p^2t^2+a_5)^\frac{a_1}{4}e^{k_pxt+\frac{a_4-2a_3}{\sqrt{4a_5-a_3^2}}\tan^{-1}\frac{2k_p^2t^2+a_3}{\sqrt{4a_5-a_3^2}}}$

$\int_{m_p}^{n_p}(k_p^4t^4+a_3k_p^2t^2+a_5)^\frac{a_1-4}{4}e^{k_pxt+\frac{a_4-2a_3}{\sqrt{4a_5-a_3^2}}\tan^{-1}\frac{2k_p^2t^2+a_3}{\sqrt{4a_5-a_3^2}}}~dt$ converges

Case $1c$: $a_5<\dfrac{a_3^2}{4}$

Then $\ln K(s)=\dfrac{a_1-4}{4}\ln(s^4+a_3s^2+a_5)+\dfrac{2a_3-a_4}{\sqrt{a_3^2-4a_5}}\tanh^{-1}\dfrac{2s^2+a_3}{\sqrt{a_3^2-4a_5}}+c_1$

$K(s)=c(s^4+a_3s^2+a_5)^\frac{a_1-4}{4}e^{\frac{2a_3-a_4}{\sqrt{a_3^2-4a_5}}\tanh^{-1}\frac{2s^2+a_3}{\sqrt{a_3^2-4a_5}}}$

$\therefore y=\int_m^nc(s^4+a_3s^2+a_5)^\frac{a_1-4}{4}e^{xs+\frac{2a_3-a_4}{\sqrt{a_3^2-4a_5}}\tanh^{-1}\frac{2s^2+a_3}{\sqrt{a_3^2-4a_5}}}~ds$

But since the above procedure in fact suitable for any complex number $s$ ,

$\therefore y_p=\int_{m_p}^{n_p}c_p((k_pt)^4+a_3(k_pt)^2+a_5)^\frac{a_1-4}{4}e^{xk_pt+\frac{2a_3-a_4}{\sqrt{a_3^2-4a_5}}\tanh^{-1}\frac{2(k_pt)^2+a_3}{\sqrt{a_3^2-4a_5}}}~d(k_pt)=k_pc_p\int_{m_p}^{n_p}(k_p^4t^4+a_3k_p^2t^2+a_5)^\frac{a_1-4}{4}e^{k_pxt+\frac{2a_3-a_4}{\sqrt{a_3^2-4a_5}}\tanh^{-1}\frac{2k_p^2t^2+a_3}{\sqrt{a_3^2-4a_5}}}~dt$

For some $x$-independent real number choices of $m_p$ and $n_p$ and $x$-independent complex number choices of $k_p$ such that:

$\lim\limits_{t\to m_p}(k_p^4t^4+a_3k_p^2t^2+a_5)^\frac{a_1}{4}e^{k_pxt+\frac{2a_3-a_4}{\sqrt{a_3^2-4a_5}}\tanh^{-1}\frac{2k_p^2t^2+a_3}{\sqrt{a_3^2-4a_5}}}=\lim\limits_{t\to n_p}(k_p^4t^4+a_3k_p^2t^2+a_5)^\frac{a_1}{4}e^{k_pxt+\frac{2a_3-a_4}{\sqrt{a_3^2-4a_5}}\tanh^{-1}\frac{2k_p^2t^2+a_3}{\sqrt{a_3^2-4a_5}}}$

$\int_{m_p}^{n_p}(k_p^4t^4+a_3k_p^2t^2+a_5)^\frac{a_1-4}{4}e^{k_pxt+\frac{2a_3-a_4}{\sqrt{a_3^2-4a_5}}\tanh^{-1}\frac{2k_p^2t^2+a_3}{\sqrt{a_3^2-4a_5}}}~dt$ converges

In fact the most difficulties for solving higher order ODEs is that it is extremely difficult to control its coefficients even applying some simple transformation types.

Try let $y=e^{kx^2}u$ ,

Then $y'=e^{kx^2}u'+2kxe^{kx^2}u$

$y''=e^{kx^2}u''+2kxe^{kx^2}u'+2kxe^{kx^2}u'+(4k^2x^2+2k)e^{kx^2}u=e^{kx^2}u''+4kxe^{kx^2}u'+(4k^2x^2+2k)e^{kx^2}u$

$y'''=e^{kx^2}u'''+2kxe^{kx^2}u''+4kxe^{kx^2}u''+(8k^2x^2+4k)e^{kx^2}u'+(4k^2x^2+2k)e^{kx^2}u'+(2kx(4k^2x^2+2k)+8k^2x)e^{kx^2}u=e^{kx^2}u'''+6kxe^{kx^2}u''+(12k^2x^2+6k)e^{kx^2}u'+(8k^3x^3+12k^2x)e^{kx^2}u$

$y''''=e^{kx^2}u''''+2kxe^{kx^2}u'''+6kxe^{kx^2}u'''+(12k^2x^2+6k)e^{kx^2}u''+(12k^2x^2+6k)e^{kx^2}u''+(24k^3x^3+36k^2x)e^{kx^2}u'+(8k^3x^3+12k^2x)e^{kx^2}u'+(16k^4x^4+48k^3x^2+12k^2)e^{kx^2}u=e^{kx^2}u''''+8kxe^{kx^2}u'''+(12k^2x^2+6k)e^{kx^2}u''+(32k^3x^3+48k^2x)e^{kx^2}u'+(16k^4x^4+48k^3x^2+12k^2)e^{kx^2}u$

$\therefore x^2(e^{kx^2}u''''+8kxe^{kx^2}u'''+(12k^2x^2+6k)e^{kx^2}u''+(32k^3x^3+48k^2x)e^{kx^2}u'+(16k^4x^4+48k^3x^2+12k^2)e^{kx^2}u)+a_1x(e^{kx^2}u'''+6kxe^{kx^2}u''+(12k^2x^2+6k)e^{kx^2}u'+(8k^3x^3+12k^2x)e^{kx^2}u)+(a_2+a_3x^2)(e^{kx^2}u''+4kxe^{kx^2}u'+(4k^2x^2+2k)e^{kx^2}u)+a_4x(e^{kx^2}u'+2kxe^{kx^2}u)+a_5x^2e^{kx^2}u=0$

$x^2(u''''+8kxu'''+(12k^2x^2+6k)u''+(32k^3x^3+48k^2x)u'+(16k^4x^4+48k^3x^2+12k^2)u)+a_1x(u'''+6kxu''+(12k^2x^2+6k)u'+(8k^3x^3+12k^2x)u)+(a_2+a_3x^2)(u''+4kxu'+(4k^2x^2+2k)u)+a_4x(u'+2kxu)+a_5x^2u=0$

$x^2u''''+(8kx^2+a_1)xu'''+(12k^2x^4+(6(a_1+1)k+a_3)x^2+a_2)u''+(32k^3x^5+(12(a_1+4)k^2+4a_3k)x^3+(2(3a_1+2a_2)k+a_4)x)u'+(16k^4x^6+4((12+2a_1)k+a_3)k^2x^4+(4(3a_1+a_2+3)k^2+2(a_3+a_4)k+a_5)x^2+2a_2k)u=0$

The issue is even much more complicated.

Answer 2

One technique for solving linear higher order ODEs is to seek solutions in the form $y(x) = y_1(x) \cdot y_2(x)$ where $y_{1,2}$ are solutions of some lower order ODE (input ODE), solutions that are known. In order to guarantee success of this method one has make sure that both the input and the target ODE (i.e. the ODE in question) have the same singular points. The target ODE has singular points at zero and at infinity and therefore we chose the Whittaker differential equation (which also has singularities at zero and infinity) as the input ODE.

Define $Q(x) := -A + k/x + (1/4-\mu^2)/x^2$. Then the equation(to be called input ODE) reads: \begin{equation} \frac{d^2 y_{1,2}(x)}{d x^2} + Q(x) y_{1,2}(x) = 0 \end{equation} is solved by \begin{eqnarray} y_{1}(x) &=& M_{\frac{\kappa}{2\sqrt{A}},\mu}(2\sqrt{A} x) \\ y_{2}(x) &=& W_{\frac{\kappa}{2\sqrt{A}},\mu}(2\sqrt{A} x) \end{eqnarray} where $M_{\cdot,\cdot}(x)$and $W_{\cdot,\cdot}(x)$ are Whittaker functions https://en.wikipedia.org/wiki/Whittaker_function . Now we differentiate our function $y(x)$. We have: \begin{eqnarray} &&y(x) = y_1(x) y_2(x) \\ &&y^{'}(x) = y_1^{'}(x) y_2(x) + y_1(x) y_2^{'}(x)\\ &&y^{''}(x) = -2 Q(x) y_1(x) y_2(x) + 2 y_1^{'}(x) y_2^{'}(x)\\ &&y^{'''}(x)=-2 Q^{'}(x) y_1(x) y_2(x) - 4 Q(x) (y_2(x) y_1^{'}(x) +y_1(x) y_2^{'}(x)) \\ &&y^{''''}(x)= (8 Q^2(x)-2 Q^{''}(x)) y_1(x) y_2(x)-6 Q^{'}(x)(y_2(x) y_1^{'}(x) +y_1(x) y_2^{'}(x)) -8 Q(x) y_1^{'}(x) y_2^{'}(x) \end{eqnarray} Now we insert the expressions above into our target ODE and recognize that there will be only three type of terms being proportional to firstly $y_1(x) y_2(x)$, secondly to $(y_1(x) y_2^{'}(x) + y_2(x) y_1^{'}(x))$ and thirdly to $y_1^{'}(x) y_2^{'}(x)$ with the proportionality constants being polynomials of order at most four in the variable $1/x$. Now all we need to do is to choose the parameters $(A,\kappa,\mu)$of our input ODE so that the polynomials in question are all identically equal to zero. Even though the later seems to be a very strong requirement it turns out that solutions do exist for a certain choice of the parameters $\left(a_i \right)_{i=1}^5$ of the target ODE.

As a matter of fact the following is true. Let $a_2$ and $a_3$ be arbitrary real numbers. Then let : \begin{eqnarray} a_1&=& 3\\ a_4&=& a_1 \cdot a_3\\ a_5&=&0 \end{eqnarray} Then the set of functions $(y_1(x) y_2(x), y_1(x)^2,y_2(x)^2)$ where \begin{eqnarray} y_{1}(x) &=& M_{\frac{\kappa}{2\sqrt{A}},\mu}(2\sqrt{A} x) \\ y_{2}(x) &=& W_{\frac{\kappa}{2\sqrt{A}},\mu}(2\sqrt{A} x) \end{eqnarray} and $(k,\mu,A)= (0,\sqrt{1-a_2}/2,-a_3/4)$ solves the target ODE. The following Mathematica code verifies that:

In[125]:= {a1, a2, a3, a4, a5} = RandomInteger[{1, 10}, 5];
a1 = 3;
a4 = a1 a3;
a5 = 0;
{k, mu, A} = {0, Sqrt[1 - a2]/2, -a3/4};
y1[x_] = WhittakerM[k/(2 Sqrt[A]), mu, 2 Sqrt[A] x];
y2[x_] = WhittakerW[k/(2 Sqrt[A]), mu, 2 Sqrt[A] x];
FullSimplify[(D[#, {x, 4}] + 
     a1/x D[#, {x, 3}] + (a3 + a2/x^2) D[#, {x, 2}] + 
     a4/x D[#, {x, 1}] + a5 #) & /@ {y1[x] y2[x], y1[x]^2, y2[x]^2}]

Out[132]= {0, 0, 0}

Update 0:

It is more likely for us to obtain new solutions if we introduce additional parameters into the system. This can be done by replacing $y_{1,2}(x) \rightarrow r_0 y_{1,2}(x) + r_1 y_{1,2}^{'}(x)$ where $r_0$ and $r_1$ are new parameters that come into play. Having done that we just repeat the procedure above and interestingly enough we find a new solution. Let $a_3$ and $a_5$ be arbitrary real numbers. Now let : \begin{eqnarray} a_1&=&3\\ a_2&=&0\\ a_4&=& \frac{3}{2} \left(a_3+\sqrt{a_3^2-4 a_5}\right) \end{eqnarray} Then the set of functions: \begin{eqnarray} \left( \begin{array}{r} (r_0 y_1(x) + r_1 y_1^{'}(x))\cdot (r_0 y_2(x) + r_1 y_2{'}(x))\\ (r_0 y_1(x) + r_1 y_1^{'}(x))\cdot (r_0 y_1(x) + r_1 y_1{'}(x))\\ (r_0 y_2(x) + r_1 y_2^{'}(x))\cdot (r_0 y_2(x) + r_1 y_2{'}(x)) \end{array} \right) \end{eqnarray} where \begin{eqnarray} y_{1}(x) &=& M_{\frac{\kappa}{2\sqrt{A}},\mu}(2\sqrt{A} x) \\ y_{2}(x) &=& W_{\frac{\kappa}{2\sqrt{A}},\mu}(2\sqrt{A} x) \end{eqnarray} and $(k,\mu,A)= (0,1/2,1/8(-a_3-\sqrt{a_3^2-4 a_5}))$ and \begin{equation} r_0=\frac{\imath \sqrt{a_5} r_1}{\sqrt{2} \sqrt{a_3-\sqrt{a_3^2-4 a_5}}} \end{equation} solves the target ODE. Again, we used Mathematica to verify our result.

In[292]:= a3 =.; a5 =.; r1 =.; r0 =.; Clear[y1]; Clear[y2];
{a1, a2, a4} = {3, 0, 3/2 (a3 + Sqrt[a3^2 - 4 a5])};
r0 = (I Sqrt[a5] r1)/(Sqrt[2] Sqrt[a3 - Sqrt[a3^2 - 4 a5]]);
{A, mu, k} = {1/8 (-a3 - Sqrt[a3^2 - 4 a5]), 1/2, 0};
y1[x_] = WhittakerM[k/(2 Sqrt[A]), mu, 2 Sqrt[A] x];
y2[x_] = WhittakerW[k/(2 Sqrt[A]), mu, 2 Sqrt[A] x];

FullSimplify[(D[#, {x, 4}] + 
     a1/x D[#, {x, 3}] + (a3 + a2/x^2) D[#, {x, 2}] + a4/x D[#, x] + 
     a5 #) & /@ {(r0 y1[x] + r1 D[y1[x], x]) (r0 y2[x] + 
      r1 D[y2[x], x]), (r0 y1[x] + r1 D[y1[x], x]) (r0 y1[x] + 
      r1 D[y1[x], x]), (r0 y2[x] + r1 D[y2[x], x]) (r0 y2[x] + 
      r1 D[y2[x], x])}]

Out[298]= {0, 0, 0}

Update 1:

Here is a third possible way to tackling the problem in question. We know that we can always annihilate the coefficient at the third order derivative by writing $y(x)=\exp(-1/4 \int(a_1/x dx)) \cdot v(x)$ and then apply the procedure described above to the equation satisfied by the function $v(x)$. If we do this we obtain another solutions to the ODE is shown below.

Solution 1:

Let $a_3$ be real and then let the following: \begin{eqnarray} a_1&=&4\\ a_2&=&0\\ a_4&=&\frac{a1 a3}{2}\\ a_5&=&0 \end{eqnarray} Then define \begin{eqnarray} y_{1}(x) &=& M_{\frac{\kappa}{2\sqrt{A}},\mu}(2\sqrt{A} x) \\ y_{2}(x) &=& W_{\frac{\kappa}{2\sqrt{A}},\mu}(2\sqrt{A} x) \end{eqnarray} where $(k,\mu,A)= (0,1/2,-a_3/4))$.

Then the set of functions $(y_1(x) y_2(x)/x, y_1(x)^2/x,y_2(x)^2/x)$ satisfy the ODE in question. The Mathematica code snippet verifies it:

In[13]:= a3 =.; a1 =.; k =.; A =.; mu =.; Clear[y1]; Clear[y2];
{a1, a2, a4, a5} = {4, 0, (a1 a3)/2, 0};
{k, A, mu} = {0, -a3/4, 1/2};
y1[x_] = WhittakerM[k/(2 Sqrt[A]), mu, 2 Sqrt[A] x];
y2[x_] = WhittakerW[k/(2 Sqrt[A]), mu, 2 Sqrt[A] x];

FullSimplify[(D[#, {x, 4}] + 
     a1/x D[#, {x, 3}] + (a3 + a2/x^2) D[#, {x, 2}] + a4/x D[#, x] + 
     a5 #) & /@ {1/x y1[x] y2[x], 1/x y1[x]^2, 1/x y2[x]^2}]

Out[18]= {0, 0, 0}    

Solution 2:

Let $a_5$ be real and then let the following: \begin{eqnarray} a_1&=&8\\ a_2&=&12\\ a_3&=&0\\ a_4&=&0 \end{eqnarray} and \begin{equation} r_0=\frac{\sqrt{a_5}}{2(-a_5)^{1/4}} r_1 \end{equation} Then define \begin{eqnarray} y_{1}(x) &=& M_{\frac{\kappa}{2\sqrt{A}},\mu}(2\sqrt{A} x) \\ y_{2}(x) &=& W_{\frac{\kappa}{2\sqrt{A}},\mu}(2\sqrt{A} x) \end{eqnarray} where $(k,\mu,A)= (0,1/2,-1/4 \sqrt{-a_5})$.

Then the set of functions: \begin{eqnarray} \left( \begin{array}{r} (r_0 y_1(x) + r_1 y_1^{'}(x))\cdot (r_0 y_2(x) + r_1 y_2{'}(x))\cdot x^{-2}\\ (r_0 y_1(x) + r_1 y_1^{'}(x))\cdot (r_0 y_1(x) + r_1 y_1{'}(x))\cdot x^{-2}\\ (r_0 y_2(x) + r_1 y_2^{'}(x))\cdot (r_0 y_2(x) + r_1 y_2{'}(x))\cdot x^{-2} \end{array} \right) \end{eqnarray} satisfy the ODE in question. The Mathematica code snippet verifies it:

In[1]:= a1 =.; a2 =.; a3 =.; a4 =.; a5 =.; r1 =.; r0 =.; Clear[y1]; \
Clear[y2];
{a1, a2, a3, a4} = {8, 12, 0, 0};
r0 = Sqrt[a5]/(2 (-a5)^(1/4)) r1;
mu = 1/2; k = 0; A = 1/4 (-Sqrt[- a5]);

y1[x_] = WhittakerM[k/(2 Sqrt[A]), mu, 2 Sqrt[A] x];
y2[x_] = WhittakerW[k/(2 Sqrt[A]), mu, 2 Sqrt[A] x];

FullSimplify[(D[#, {x, 4}] + 
     a1/x D[#, {x, 3}] + (a3 + a2/x^2) D[#, {x, 2}] + a4/x D[#, x] + 
     a5 #) & /@ {(r0 y1[x] + 
      r1 D[y1[x], x]) (r0 y2[x] + r1 D[y2[x], x])/
     x^2, (r0 y1[x] + r1 D[y1[x], x]) (r0 y1[x] + r1 D[y1[x], x])/
     x^2, (r0 y2[x] + r1 D[y2[x], x]) (r0 y2[x] + r1 D[y2[x], x])/x^2}]

Out[7]= {0, 0, 0}

Of course what we have found are only very particular solutions to the target ODE however one can readily see that by introducing additional transformations and in turn additional parameters we might be able to find more solutions. One possibility to try would be to firstly try a change of variables $x\rightarrow f(x)$ and $d/d x \rightarrow 1/f^{'}(x) d/d x$ in the input equation and then choose $f(x)$ in such a way that the target equation is solved.

Answer 3

A power series method might do the trick. Substitute $$ y = \sum_{k=0}^\infty c_k x^k$$ in the differential equation $$x^2 y^{(4)} + a_1 x y^{(3)} + \cdots + a_5 x^2 y = 0.$$ You will obtain a series in which the coefficients are difference equations in the unknowns $c_k$. For concrete constants $a_n$, this difference equation maybe can be solved.

Answer 4

Since our ODE is of fourth order and has regular singularities at zero and infinity it would make sense the map it onto a different ODE with the same properties and with known solutions. In here the generalized hypergeometric function $f(x):=F_{0,3}[b_1,b_2,b_3;x]$ is a good candidate. That function satisfies the following ODE (see https://en.wikipedia.org/wiki/Generalized_hypergeometric_function): \begin{equation} f(x)= \frac{d}{d x} \prod\limits_{n=1}^3 \left( x \frac{d}{d x} + b_n-1 \right) f(x) \end{equation} Now let us define a scaling function $\theta(x) := A x^r$ and in the ODE above let us change the variable as follows: $x \rightarrow \theta(x)$ and $d/d x \rightarrow 1/\theta^{'}(x) d/d x$. Then it turns out that the rescaled function satisfies the following ODE: \begin{eqnarray} 0=-A r^4 f(x) x^{r-4}+\frac{((b_1+b_2+b_3-3) r+6) f^{(3)}(x)}{x}+\frac{\left((b_2 (b_3-2)-2 b_3+b_1 (b_2+b_3-2)+3) r^2+3 (b_1+b_2+b_3-3) r+7\right) f''(x)}{x^2}+\frac{((b_1-1) r+1) ((b_2-1) r+1) ((b_3-1) r+1) f'(x)}{x^3}+f^{(4)}(x) \end{eqnarray} where with a slight abuse of notation we can write $f(x):=f(\theta(x))$. By comparing the above with our original ODE we write: \begin{eqnarray} r&=&4\\ -A r^4 &=& a_5\\ (b_1-1) r+1 &=& 0\\ \left((b_2 (b_3-2)-2 b_3+b_1 (b_2+b_3-2)+3) r^2+3 (b_1+b_2+b_3-3) r+7\right)&=& a_2\\ ((b_1+b_2+b_3-3) r+6) &=& a_1 \end{eqnarray} This has a following solution: \begin{eqnarray} r&=&4\\ A&=&-\frac{a_5}{256}\\ b_1 &=& \frac{3}{4}\\ b_2 &=& \frac{1}{8} \left( 3+a_1- \sqrt{(1-a_1)^2-4 a_2}\right)\\ b_3 &=& \frac{1}{8} \left( 3+a_1+ \sqrt{(1-a_1)^2-4 a_2}\right) \end{eqnarray} and therefore a function \begin{equation} v(x):= F_{0,3}\left[b_1,b_2,b_2,-\frac{a_5}{256} x^4\right] \end{equation} solves the equation: \begin{equation} 0=\frac{a_1 v^{(3)}(x)}{x}+\frac{a_2 v''(x)}{x^2}+a_5 v(x)+v^{(4)}(x) \end{equation} which is the original ODE subject to $a_3=$ and $a_4=0$.

The code snippet verifies the result:

In[53]:= 
A1 =.; A2 =.; A3 =.; A4 =.; A5 =.;
{A, r} = {-A5/256, 4};
{b[1], b[2], b[3]} = {3/4, 1/8 (3 + A1 - Sqrt[(1 - A1)^2 - 4 A2]), 
   1/8 (3 + A1 + Sqrt[(1 - A1)^2 - 4 A2])};
v[x_] = HypergeometricPFQ[{}, {b[1], b[2], b[3]}, A x^r];
FullSimplify[
 A5 v[x] + (A2/x^2) D[v[x], {x, 2}] + A1 /x D[v[x], {x, 3}] + 
  D[v[x], {x, 4}]]

Out[57]= 0

Update: Since mapping generalized hypergeometric functions onto our ODE appears to be quite fruitful let us try now with another specimen of this kind. We take $f(x):=F_{2,3}\left[a_1,a_2;b_1,b_2,b_3;x\right]$. We take $(p,q)=(2,3)$ and our function satisfies the following ODE: \begin{equation} x \prod\limits_{n=1}^p \left(x \frac{d}{d x} + a_n\right) f(x)= x \frac{d}{d x} \prod\limits_{n=1}^q \left( x \frac{d}{d x} + b_n-1 \right) f(x) \end{equation} and now we repeat the procedure above. This means we firstly rescale the abscissa $x \rightarrow A x^r$ and then rescale the ordinate $f(x) \rightarrow x^s v(x)$. As usual we carry out those calculations in Mathematica because the formulae are unwieldy otherwise. In this case we have eight parameters to tweek firstly the five parameters of the hypergeometric function $(a_j)_{j=1}^p$, $(b_j)_{j=1}^q$ and the two parameters of the asbcissa scaling $A$,$r$ and finally the ordinate scaling exponent $s$. It turns out that the following choice is fruitful. We have:

\begin{eqnarray} s&=&\frac{1}{2}\left(-5+A_1 + \sqrt{(1-A_1)^2-4 A_2} \right)\\ a_1&=& -\frac{s}{2}\\ a_2&=&\frac{A_4-A_3(1+s)}{2 A_3}\\ b_1&=&\frac{-s+2}{2}\\ b_2&=&\frac{1-s}{2}\\ b_3&=&\frac{A_1-3-2 s}{2} \end{eqnarray} and \begin{eqnarray} A&=&-\frac{A_3}{4}\\ r&=&2 \end{eqnarray} and therefore the function \begin{eqnarray} v(x)&=& x^{-s} F_{2,3}\left[a_1,a_2;b_1,b_2,b_3;A x^r\right] \end{eqnarray} solves the equation: \begin{equation} 0=\frac{A_1 v^{(3)}(x)}{x}+\left(\frac{A_2}{x^2}+A_3\right) v''(x)+\frac{A_4 v'(x)}{x}+v^{(4)}(x) \end{equation} which is the original equation subject to the coefficient at the zeroth derivative (the right-most term) being equal to zero.

Again the code snippet provides a verification:

In[61]:= A1 =.; A2 =.; A3 =.; A4 =.; A5 =.;
s = 1/2 (-5 + A1 + Sqrt[(-1 + A1)^2 - 4 A2]);
{a[1], a[2], b[1], b[2], b[3]} = {-s/2, (A4 - A3 (1 + s))/(
   2 A3), (-s + 2)/2, (1 - s)/2, (A1 - 3 - 2 s)/2};
{A, r} = {-A3/4, 2};
v[x_] = x^(-s) HypergeometricPFQ[{a[1], a[2]}, {b[1], b[2], b[3]}, 
    A x^r];
FullSimplify[
 A4 /x D[v[x], x] + (A3 + A2/x^2) D[v[x], {x, 2}] + 
  A1 /x D[v[x], {x, 3}] + D[v[x], {x, 4}]]


Out[66]= 0

Answer 5

This is a follow-up to the answer given by doraemonpaul above.Let us express the unknown function through its Laplace transform as follows: \begin{equation} y(x)=\int\limits_0^\infty e^{-x s} K(s) ds \end{equation} By inserting the above into the ODE and then by integrating by parts in order to eleiminate the powers of $x$ from the integrand in the $s$-domain we obtain the following equation: \begin{eqnarray} &&\frac{d^2}{d s^2} \left[s^4 K(s) \right] - a_1 \frac{d}{d s} \left[s^3 K(s)\right]+(a_2+a_3 \frac{d^2}{d s^2}) \left[s^2 K(s) \right] - a_4 \frac{d}{d s}\left[ s K(s) \right] + a_5 \frac{d^2}{d s^2} \left[ K(s) \right] =\\ && -a_5 x K(0) - a_5 K^{'}(0) \end{eqnarray}

Now let us for simplicity assume that $K(0)=K^{'}(0)=0$ which means that $y(x) = O(1/x^3)$ when $x\rightarrow \infty$. Then we have: \begin{eqnarray} &&\frac{a_2 s^2 K(s) + \frac{d}{d s}\left( \frac{d}{d s}\left((a_5+a_3 s^2+s^4)K(s)\right) - (a_1 s^3 + a_4 s)K(s)\right)}{a_5+a_3 s^2+s^4}=\\ && \frac{K[s] \left(s^2 (-3 a_1+a_2+12)+2 a_3-a_4\right)}{a_3 s^2+a_5+s^4}-\frac{s K'(s) \left((a_1-8) s^2-4 a_3+a_4\right)}{a_3 s^2+a_5+s^4}+K''(s) =0 \end{eqnarray}

Then by substituting for $u=s^2$, $d/d s = 2 \sqrt{u} d/d u$ and $d^2/d u^2 = 2 d/d u + 4 u d^2/d u^2$ we get a following ODE satisfied by the Laplace transform: \begin{equation} \frac{K[u] (u (-3 a_1+a_2+12)+2 a_3-a_4)}{4 u (u (a_3+u)+a_5)}+\frac{K'(u) (a_5-u ((a_1-9) u-5 a_3+a_4))}{2 u (u (a_3+u)+a_5)}+K''(u)=0 \end{equation} Now we use the common technique to eliminate the coefficient at the first derivative. We write $K(u) = M(u) \cdot V(u)$ where $M(u)$ is the exponential from minus one half times the antiderivative of the coefficient in question. In this case it reads: \begin{equation} M(u):=\frac{(u (a_3+u)+a_5)^{\frac{a_1-8}{8}} \exp \left(-\frac{(a_1 a_3-2 a_4) \tan ^{-1}\left(\frac{a_3+2 u}{\sqrt{4 a_5-a_3^2}}\right)}{4 \sqrt{4 a_5-a_3^2}}\right)}{\sqrt[4]{u}} \end{equation} Now define: \begin{eqnarray} P^{(4)}(u)&:=&u^4 \left(-a_1^2+2 a_1+4 a_2+3\right)+\\ &&u^3 (a_3 (-2 a_1+4 a_2+6)-2 (a_1-3)a_4)+\\ &&u^2 \left(2 a_5 (-3 a_1+2 a_2+3)+3 a_3^2+2 a_3 a_4-a_4^2\right)+\\ &&u (6 a_3 a_5-2 a_4 a_5)+\\ &&3 a_5^2 \end{eqnarray} and then the function $V(u)$ satisfies the following ODE: \begin{equation} V^{''}(u)+\frac{P^{(4)}(u)}{16 u^2(a_5+a_3 u + u^2)^2}V(u)=0 \end{equation} which maps to the Heun equation https://en.wikipedia.org/wiki/Heun_function .