# Preface

This secton is devoted to a special class of variable coefficients linear differential equations. This remarkable class was discovered by Leonhard Euler who showed that these differential equations could be solved explicitely via elementary functions.

# Euler Equations

An Euler equation (also known as the Euler-Cauchy equation, or equidimensional equation) is a linear homogeneous ordinary differential equation with variable coefficients of the following form:

$a_n x^n y^{(n)} + a_{n-1} x^{n-1} y^{(n-1)} + \cdots + a_1 x\, y' + a_0 y = f(x) ,$
where the coefficients $$a_0 , a_1 , \ldots , a_n$$ are constants and the driving function f(x) is given. Since the leading term in the Euler equation is a multiple of xn, the origin should be excluded from consideration. The Euler equation is usually considered on semi-infinite intervals either $$(0, \infty ) \quad\mbox{or}\quad (-\infty , 0) .$$ Let
$L\left[ x, \texttt{D} \right] = a_n x^n \texttt{D}^n + a_{n-1} x^{n-1} \texttt{D}^{n-1} + \cdots + a_1 x\,\texttt{D} + a_0 , \qquad \texttt{D} = \frac{{\text d}}{{\text d}x} ,$
be the linear differential operator corresponding to the Euler equation. We start with finding functions annihilated by the Euler operator; in other words, we are going to determine the fundamental set of solutions for the homogeneous equation $$L\left[ x, \texttt{D} \right] y = 0$$ on the positive semi-axis $$(0, \infty ) .$$

There are known two methods to solve the homogeneous Euler equation. One is to try a power function $$y = x^m ,$$ which leads to
$\frac{{\text d}^k}{{\text d}x^k} \, x^m = m^{\underline{k}} \, x^{m-k} = m(m-1) \cdots (m-k+1)\, x^{m-k} , \quad k =1,2, \ldots , n .$
Here $$\displaystyle m^{\underline{k}} = m(m-1) \cdots (m-k+1)$$ is k-th falling factorial of number m. Applying the above Euler operator to a trial solution xm, we obtain
$L\left[ x, \texttt{D} \right] x^m = x^m \left( a_n m^{\underline{n}} + a_{n-1} m^{\underline{n-1}} + \cdots + a_1 m + a_0 \right) .$
Canceling common multiple xm yields an auxiliary algebraic equation
$P(m) = a_n m^{\underline{n}} + a_{n-1} m^{\underline{n-1}} + \cdots + a_1 m + a_0 =0 ,$
where P(m) is an auxiliary polynomial of degree n (in accordance to the degree of the Euler operator). If m is a root of the above algebraic equation, then $$y = x^m$$ is a solution of the n-th order Euler homogeneous equation. We postpone analyzing the fundamental set of solutions, which depends on whether the roots of the auxiliary algebraic equation are real or complex and whether some of them have multiplicities others than 1.

Because of the particularly simple equidimensional structure, the Euler equation can be replaced with an equivalent constant coefficients differential equation, which can then be solved explicitly. Using change of independent variable $$x = e^t , \quad t = \ln x ,$$ it follows that

\begin{align*} \frac{{\text d} y}{{\text d}x} &= \frac{{\text d} y}{{\text d}t} \, \frac{{\text d} t}{{\text d}x} = \frac{1}{x} \, \frac{{\text d} y}{{\text d}t} , \\ \frac{{\text d}^2 y}{{\text d}x^2} &= \frac{1}{x} \, \frac{{\text d}}{{\text d}x} \left( \frac{{\text d} y}{{\text d}t} \right) + \frac{{\text d} y}{{\text d}t} \left( - \frac{1}{x^2} \right) \qquad \mbox{product rule and chain rule} \\ &= \frac{1}{x} \left( \frac{1}{x}\, \frac{{\text d}^2 y}{{\text d}t^2} \right) + \frac{{\text d} y}{{\text d}t} \left( - \frac{1}{x^2} \right) = \frac{1}{x^2} \left( \frac{{\text d}^2 y}{{\text d}t^2} - \frac{{\text d} y}{{\text d}t} \right) , \end{align*}
and so on. So we see that every derivative with respect to x, upon transformation to a new variable t, has a multiple reciprocal to x; therefore, all these power multiples are canceled out and we obtain a constant coefficient differential equation. We summarize our observation in the table:
Solutions of P(m) = 0  Fundamental set of solutions
Simple roots m1, m2  $$y_1 = x^{m_1} , \ y_2 = x^{m_2} , \ldots$$
Double roots $$m_1 = m_2$$  $$y_1 = x^{m_1} , \ y_2 = x^{m_1} \ln x$$
Triple roots $$m_1 = m_2 = m_3$$  $$y_1 = x^{m_1} , \ y_2 = x^{m_1} \ln x , \ y_3 = x^{m_1} \left( \ln x \right)^2$$
Complex roots $$m_1 = \overline{m_2} = \alpha + \beta {\bf j}$$  $$y_1 = x^{m_1} \cos \left( \beta\,\ln x \right) , \ y_2 = x^{m_1} \sin \left( \beta\,\ln x \right)$$

To make our exposition transparent, we knock out second order Euler equations and then consider the third order equations.

Second Order Euler Equations

Consider second order homogeneous Euler equation

$a\,x^2 y'' + b\,x\, y' + c\,y =0 ,$
where a, b, and c are some given real numbers. For a trial solution $$y= x^m ,$$ we have
\begin{align*} y' &= \frac{{\text d}y}{{\text d}x} = m\, x^{m-1} , \\ y'' &= \frac{{\text d}^2 y}{{\text d}x^2} = m(m-1)\, x^{m-2} . \end{align*}
Plugging into the Euler equation, we get
\begin{align*} 0 &= a\,x^2 y'' + b\,x\,y' + c\,y \\ &= a\, x^2 m(m-1) \,x^{m-2} + b\,x\,m\,x^{m-1} + c\,x^m \\ &= a\, m(m-1)\,x^m + b\,m\,x^m + c\,x^m \\ &= \left[ a\,m(m-1) + b\,m + c \right] x^m . \end{align*}
Thus, $$y= x^m$$ is a solution of the Euler equation if and only if m is a root of the auxiliary polynomial
$P(m) = a\,m(m-1) + b\,m + c ,$
which simplifies to
$P(m) = a\,m^2 + \left( b-a \right) m + c .$
If P(m) has two distinct roots, either real or complex, we know how to solve the equation. To recall briefly, if P(m) has two distinct real roots m1 and m2, then $$x^{m_1} \quad\mbox{and}\quad x^{m_2}$$ are two linearly independent solutions of the Euler equation; hence, the general solution becomes
$y= C_1 x^{m_1} + C_2 x^{m_2} ,$
for some arbitrary constants C1 and C2.

If P(m) has two distinct complex roots $$m_1 \quad\mbox{and}\quad m_2 ,$$ they must be complex conjugates, so we can write $$m_1 = \alpha + \beta {\bf j}$$ and $$m_2 = \alpha - \beta {\bf j} ,$$ where α and β are real numbers and j is a unit vector in positive vertical direction on the complex plane, so $${\bf j}^2 =-1 .$$

It's perhaps not so clear what a complex power of x means, but since x is assumed to be positive, we have $$x = e^{\ln x} = e^t .$$ If the complex power makes any kind of sense, we ought to have

\begin{align*} x^{\alpha + \beta \,{\bf j}} &= \left( e^{\ln x} \right)^{\alpha + {\bf j}\,\beta} = e^{\ln x \left( \alpha + {\bf j}\,\beta \right)} \\ &= e^{\alpha\,\ln x + {\bf j}\,\beta\,\ln x} = e^{\alpha\,\ln x} \cdot e^{{\bf j}\,\beta\,\ln x} \\ &= \left( e^{\ln x} \right)^{\alpha} \, e^{{\bf j}\,\beta\,\ln x} \qquad\mbox{Euler formula} \\ &= x^{\alpha} \left[ \cos \left( \beta\,\ln x \right) + {\bf j}\,\sin \left( \beta\,\ln x \right) \right] . \end{align*}
Thus, we define
$x^{\alpha + {\bf j}\,\beta} = x^{\alpha} \, \cos \left( \beta\,\ln x \right) + {\bf j}\,x^{\alpha} \,\sin \left( \beta\,\ln x \right) .$
To justify this definition, we need to check
$\frac{{\text d}}{{\text d}x} \, x^{\alpha + {\bf j}\,\beta} = \left( \alpha + {\bf j}\,\beta \right) x^{\alpha -1 + {\bf j}\,\beta} \qquad\mbox{and} \qquad \overline{x^{\alpha + {\bf j}\,\beta}} = x^{\alpha - {\bf j}\,\beta}$
Verifications of these equations we leave to Mathematica.
f[x_, a_] = x^a *(Cos[b*Log[x]] + I*Sin[b*Log[x]])
TrueQ[Simplify[D[f[x, a], x]] - (a + b*I)*f[x, a - 1] == 0]
True
We can also check that our two complex solutions $$y_1 = x^{\alpha + {\bf j}\,\beta}$$ and $$y_2 = x^{\alpha - {\bf j}\,\beta}$$ are linearly independent by computing their Wronskian
$W[y_1 , y_2 ] = \det \begin{bmatrix} y_1 & y_2 \\ y'_1 & y'_2 \end{bmatrix} = -2{\bf j}\,\beta \,x^{2\alpha -1} \ne 0 \quad \mbox{for}\quad x > 0 .$
g[x_] = Refine[
x^a *(Cos[b*Log[x]] + I*Sin[b*Log[x]]), {Element[a, Reals], Element[b, Reals], x > 0}]
g1[x_] = Refine[
Conjugate[g[x]], {Element[a, Reals], Element[b, Reals], x > 0}]
Wronskian[{g[x], g1[x]}, x]
-2 I b x^(-1 + 2 a)
Therefore, we conclude that in the case of complex roots the general complex solution of the Euler equation is
$y(x) = A_1 x^{\alpha + {\bf j}\,\beta} + A_2 x^{\alpha - {\bf j}\,\beta} ,$
where A1 and A2 are arbitrary complex numbers. Each complex number is an ordered pair of two real numbers, so the above complex solution contains four arbitrary real constants. Since we are after a real-valued solution, we have to exclude two of them and impose additional constraint: $$A_2 = \overline{A_1} ,$$ so they should be complex conjugate. In this case, we reduce our general form of solution to the form
$y(x) = C_1 x^{\alpha} \, \cos \left( \beta \, \ln x \right) + C_2 x^{\alpha} \, \sin \left( \beta \, \ln x \right) ,$
for some real constants C1 and C2.

Finally, we consider the case when the auxiliary polynomial P(m) has one double root: $$P(m) = a\left( m- m_1 \right)^2 .$$ Then we have the only one solution expressed through a power function: $$y_1 = x^{m_1} .$$ So we we need to find a second solution that is independent of y1. This can be done by the method of reduction of order. To apply this method we look for a solution of the Bernoulli form $$y = u(x)\,y_1 ,$$ where y1 is the solution we already know. Before plugging this in, we need to rewrite the differential equation in the right form by expanding the auxiliary polynomial: $$P(m) = a \left( m- m_1 \right)^2 = a\, m^2 - 2a\, m_1 m + a\,m_1^2 .$$ Equating it to $$P(m) = a\,m^2 + \left( b-a \right) m + c ,$$ we find

$-2a\, m_1 = b-a, \qquad c= a\, m_1^2 \qquad\Longrightarrow \qquad b = a \left( 1 - 2 m_1 \right) .$

Finally, there is the case where P(m) has only one real root r of multiplicity two. In this case, we get one solution $$y_1 = x^r$$ of the Euler equation and we need to find a second solution that is independent of y1. This can be done by the method of reduction of order. To apply this method we look for a solution of the form $$y = u(x)\,y_1 = u(x)\,x^r .$$ Before plugging this in, we need to rewrite the differential equation in the right form.

We need to express coefficients in terms of the double root:

$P(m) = a\left( m-r \right)^2 = a\,m^2 + \left( b-a \right) m + c = a\, m^2 -2a\,m\,r + a\,r^2 .$
Thus, we must have
\begin{align*} c &= a\, r^2 , \\ -2a\,r = b-1 \qquad &\Longrightarrow \qquad b = \left( 1 - 2r \right) . \end{align*}
Plugging these values into the Euler equation, we can rewrite it as
$a\, x^2 y'' + a\,x \left( 1 - 2r \right) y' + cr^2 y =0 \qquad \Longrightarrow \qquad x^2 y'' + x\left( 1 - 2r \right) y' + r^2 y =0 .$
Differentiating our trial solution u xr, we obtain
\begin{align*} y' &= u'\, x^r + r\,u\,x^{r-1} , \\ y'' &= u'' \, x^r + 2r\,u\,x^{r-1} + r \left( r-1 \right) u\, x^{r-2} . \end{align*}
Now plug these derivatives into equation $$x^2 y'' + x\left( 1 - 2r \right) y' + r^2 y =0$$ to obtain
\begin{align*} 0 &= x^2 y'' + \left( 1 - 2r \right) y' + r^2 y \\ &= x^2 \left[ u'' \, x^r + 2r\,u\,x^{r-1} + r \left( r-1 \right) u\, x^{r-2} \right] + x \left( 1 - 2r \right) \left[ u'\, x^r + r\,u\,x^{r-1} \right] + r^2 u\, x^r \\ &= u'' x^{r+2} + 2r\,u' x^{r+1} + r \left( r-1 \right) u\,x^r + \left( 1 - 2r \right) u' x^{r+1} + r \left( 1 - 2r \right) u\, x^r + u\,r^2 x^r \\ &= u'' x^{r+2} + u' \left[ 2r + \left( r-1 \right) \right] x^{r+1} + u \left[ r \left( r-1 \right) + r \left( 1 - 2r \right) + r^2 \right] x^r \\ &= u'' x^{r+2} + u' x^{r+1} . \end{align*}
Thus, we get the equation where dependent variable is missing:
$u'' x^{r+2} + u' x^{r+1} =0 .$
If we set $$v= u' ,$$ we get a separable equation for v:
$v' x^{r+2} + v\, x^{r+1} =0 \qquad\Longrightarrow \qquad \frac{v'}{v} = -\frac{1}{x} .$
Integration yields
$\frac{{\text d}v}{v} = -\frac{{\text d}x}{x} \qquad\Longrightarrow \qquad \ln v = -\ln x + C = \ln x^{-1} +C = \ln C\,x^{-1} ,$
where C is a constant of integration. Therefore, we have
$v = C\, x^{-1} .$
We’re only looking for one u that works, so we can choose the value of C, say C = 1. This gives us
$u' = v = x^{-1} \qquad\Longrightarrow \qquad u = \ln x + C .$
Again, we can set the constant of integration to 0 because we’re only looking for one u that works. Plugging into our trial solution, $$y = u(x)\,y_1 = u(x)\,x^r ,$$ we conclude that the second solution of the Euler equation in the case of a double root r is $$y = x^r \ln x , \quad x>0.$$ It’s pretty clear that these are independent solutions, so the general solution of the Euler equation in the case of a double root r is
$y = v = C_1 x^r + C_2 x^r \ln x .$

Example: Consider the homogeneous differential equation

$x^2 y'' + 7x \,y' -7\,y =0 .$
Guided by the theory, we only need to find two linearly independent solutions of the above equation. Substituting the trial solution $$y = x^m ,$$ we obtain auxiliary algebraic equation
$m^2 + 6\, m - 7 =0 .$
Since it has two distinct real roots $$m = -7 \quad\mbox{and} \quad m =1 ,$$ we get the general solution
$y = C_1 x^{-7} + C_2 x .$

Example: Solve the homogeneous Euler equation

$4x^2 y'' +8x \,y' + y =0 .$
Upon substitution of the trial solution $$y = x^m ,$$ we obtain auxiliary algebraic equation
$4\,m^2 + 4\, m + 1 =0 \qquad \mbox{or}\qquad \left( 2m+1 \right)^2 =0,$
which has one double root m = -1/2. The general solution becomes

$y = C_1 x^{-1/2} + C_2 x^{-1/2} \, \ln x .$

Example: Solve the initial value problem

$x^2 y'' -x \,y' + 5\, y =0 , \qquad y(1) =1, \quad y' (1) =5.$
The substitution $$y = x^m ,$$ yields
$x^2 m \left( m-1 \right) x^{m-2} -x \,m \, x^{m-1} + 5\, x^m = x^m \left( m^2 -2m +5 \right) = 0 when $$m^2 -2m +5 =0 .$$ From the quadratic formula we find that the roots are $$m_1 = 1 + 2{\bf j}$$ and $$m_2 = 1 - 2{\bf j} .$$ With the identifications α = 1 and β = 2, we find the general solution: \[ y = x \left[ C_1 \cos (2\,\ln x) + C_2 \sin (2\,\ln x) \right] .$
By applying the initial conditions $$y(1) =1, \ \ y'(1) =5$$ to the foregoing solution and using $$\ln 1 =0 ,$$ we then find, in turn, that $$C_1 =1$$ and $$C_2 =2 .$$ Hence, the solution of the initial value problem is
$y = x \left[ \cos (2\,\ln x) + 2\, \sin (2\,\ln x) \right] .$
Finally, we plot the solution
f[x_] = x*(Cos[2*Log[x]] + 2*Sin[2*Log[x]])
Plot[f[x], {x, 1, 100}, PlotStyle -> Thick,
PlotLabel -> Style["x [sin(2 ln x) + 2 sin(2 ln x)]", FontSize -> 16],
Background -> LightYellow]

Example: Solve the equation

$x^2 y'' -3x \,y' + 3\,y = 3\,\ln x \qquad (x>0) .$
With the substitution $$x = e^t \ \mbox{or} \ t = \ln x ,$$ it follows

\begin{align*} \frac{{\text d}y}{{\text d}x} &= \frac{{\text d}y}{{\text d}t} \,\frac{{\text d}t}{{\text d}x} = \frac{1}{x} \, \frac{{\text d}y}{{\text d}t} = \frac{1}{x}\, \dot{y} \\ \frac{{\text d}^2 y}{{\text d}x^2} &= \frac{1}{x^2} \left( \frac{{\text d}^2 y}{{\text d}t^2} - \frac{{\text d}y}{{\text d}t} \right) = \frac{1}{x^2} \left( \ddot{y} - \dot{y} \right) . \end{align*}
Substituting in the given differential equation and simplifying yields
$\ddot{y} -4\,\dot{y} + 3\,y = 3\,t .$
Since the last equation has constant coefficients, the characteristic polynomial corresponding to the homogeneous equation is $$\lambda^2 -4\,\lambda + 3 = \left( \lambda -3 \right) \left( \lambda -1 \right) .$$ Since its nulls are λ =1 and &lambda: =3, the general solution in variable t of the corresponding homogeneous equation is
$y_h = C_1 e^t + C_2 e^{3t} ,$
with some constants C1 and C2.

Using method of undetermined coefficients, we seek a particular solution as a polynomial of first degree (similar to the right-hand side function):

$y_p = a+ bt ,$
This assumption leads to b =1 and a = 4/3. Using $$y = y_h + y_p ,$$ we get the general solution in variable t:
$y = y_h + y_p = C_1 e^t + C_2 e^{3t} + \frac{4}{3} + t .$
$y(x) = C_1 x + C_2 x^3 + \frac{4}{3} + \ln x .$

Example: Solve the nonhomogeneous Euler equation for positive x:

$x^2 y'' -x\,y' + y = x^3 e^x .$
Since the given equation is nonhomogeneous, we first need to find the fundamental set of solutions for the associated homogeneous equation $$x^2 y'' -x\,y' + y = 0 .$$ Using a trial solution $$y = x^m ,$$ we determine the auxiliary equation $$m(m-1) -m + 1=0$$ or $$(m - 1)^2 =0 .$$ In case of one double root, we obtain two linearly independent solutions $$y_1 = x \quad\mbox{and}\quad y_2 = x\,\ln x .$$ Then the general solution of the associated homogeneous equation becomes $$y_h = C_1 y_1 + C_2 y_2 = C_1 x + C_2 x \,\ln x .$$

Now we are ready to attack our main nonhomogeneous equation. According to the variation of parameters, we seek a particular solution as a linear combination:

$y_p = A(x)\, y_1 + B(x) \, y_2 = A(x)\, x + B(x) \,x\, \ln x ,$
where derivatives of the unknown functions A(x) and B(x) are determined from the Lagrange system of algebraic equations
\begin{align*} A' (x) \, x + B' (x) \,x\,\ln x &= 0, \\ A' (x) + B' (x) \, \ln x + B' (x) &= x\, e^x . \end{align*}
Solving this system is not hard to obtain
$\begin{split} A' (x) &= -x\,e^x \, \ln x , \\ B' (x) &= x\, e^x ; \end{split} \qquad\Longrightarrow \qquad \begin{split} A (x) &= -\int x\,e^x \, \ln x \,{\text d}x = e^x - \mbox{Ei} (x) - e^x \left( x-1 \right) \ln x , \\ B' (x) &= \int x\, e^x \,{\text d}x = e^x \left( x-1 \right) , \end{split}$
where Ei(x) is the exponential integral, which has to be understood in terms of the Cauchy principal value due to the singularity of the integrand at zero:
$\mbox{Ei} (x) = - (\mbox{V.P.}) \,\int_{-x}^{\infty} \frac{e^{-t}}{t} \, {\text d}t .$
Using the expressions for functions A(x) and B(x), we obtain a particular solution
$y_p = x\,e^x -x\,\mbox{Ei} (x) .$
Finally, the general solution of the given nonhomogeneous Euler equation becomes
$y = y_h + y_p = C_1 x + C_2 x \,\ln x + x\,e^x -x\,\mbox{Ei} (x) .$
Finally we plot two solutions for C1 = 1 and C2 = 2 (blue curve) or -2 (brown curve).
c1 = 1; c2 = -2;
g2[x_] = c1*x + c2*x*Log[x] + x*Exp[x] - x*ExpIntegralEi[x]
Plot[{g1[x], g2[x]}, {x, 0.01, 1.6},
PlotStyle -> {{Thick, Blue}, {Thick, Brown}},
PlotLabel -> Style["Two solutions", FontSize -> 16],
Background -> LightYellow]

Example: Suppose that we know the fundamental set of solutions $$y_1= x^3 \quad \mbox{and} \quad y_2 =x^{-1}$$ for the differential equation $$x^2 y'' -x\,y' -3y =0 .$$ If the right-hand side function is
$$(x-3)^2 /x ,$$ we need to solve the nonhomogeneous equation

$L[x, \texttt{D} ]\, y =(x-3)^2/x , \qquad\mbox{where} \quad L[x, \texttt{D} ] = x^2 \texttt{D}^2 -x\,\texttt{D} -3 .$
Here $$\texttt{D} = {\text d}/{\text d}x$$ is the derivative operator. The second order homogeneous differential equation can be written as
$y'' +p(x)y' +q(x)y=0 .$

Then

$p(x) = - W'/W ; \qquad q(x) =(y'_1 y''_2 -y''_1 y'_2)/W ,$
where $$W(x)$$ is the Wronskian of $$y_1 , y_2 .$$

Now we determine the linear differential operator $$L[x, \texttt{D}]$$ having $$y_1 (x) \ \mbox{and} \ y_2 (x)$$ as solutions of $$L[x, \texttt{D}]\, y=0.$$

y1[x_] = x^3
y2[x_] = 1/x
w[x_] = Wronskian[{y1[x], y2[x]}, x]
Out[3]= -4x
p[x_] = -w'[x]/w[x]
Out[4]= -(1/x)
q[x_] = (y1'[x] y2''[x] - y1''[x] y2'[x])/w[x]
Out[5]= -(3/x^2)
L[x_, y_] = y''[x] + p[x]*y'[x] + q[x]*y[x] (* operator *)
Out[6]= -((3 y[x])/x^2) - Derivative[1][y][x]/x + (y^$Prime]\[Prime])[x] f[x_] = (x - 3)^2 /x; basis = {x^3, 1/x} W[x_] = Table[D[basis, {x, i}], {i, 0, 1}]// Flatten Out[9]= {x^3, 1/x}, {3 x^2, -(1/x^2)} F[x_] = {0, f[x]} Out[10]= {0, (-3 + x)^2/x} v[x_] = Integrate[Inverse[W[x]].F[x], x] Out[11]= {1/4 (-(9/(2 x^2)) + 6/x + Log[x]), 1/4 (-((9 x^2)/2) + 2 x^3 - x^4/4)} yp[x_] = Simplify[basis.v[x]] (* a particular solution *) Out[12]= 1/16 x (-36 + 32 x - x^2 + 4 x^2 Log[x]) Simplify[L[x, yp] == f[x]] (* check solution *) Out[13]= True Third Order Euler Equations For arbitrary real numbers a, b, c, and d, the auxiliary polynomial \[ P(m) = a\,m \left( m-1 \right) \left( m-2 \right) + b\,m \left( m-1 \right) + c\,m + d = a\, m^3 + m^2 \left( b - 3a \right) + m \left( c-b + 2a \right) + d \qquad (a \ne 0),$
corresponding to the third order Euler equation

$a\,x^3 y''' + b\,x^2 y'' + c\,x\,y' + d\,y =0 \qquad (a \ne 0) ,$
has always one real null. Two other nulls could be either complex conjugate pair (recall that we consider only equations with real coefficients) or two real numbers. Depending on these cases, the general solution is expressed as a linear combination through the corresponding fundamental set of solutions. However, the pattern is clear. Therefore, we conclude this section with two illustrative examples.

Example: Consider a third order Euler equation

$x^3 y''' + x\,y' - y =0 .$
Plugging a trail solution $$y = x^m$$ into the given equation, we obtain the auxiliary polynomial
$P(m) = m^3 -3\, m^2 + 3\,m -1 = \left( m-1 \right)^3 ,$
which has one triple root m = 1. Therefore, we get three linearly independent solutions
$y_1 = x, \quad y_2 = x\,\ln x , \quad y_3 = x \left( \ln x \right)^2 ,$
through which the general solution is expressed:
$y = C_1 y_1 + C_2 y_2 + C_3 y_3 = C_1 x + C_2 x\,\ln x + C_3 x \left( \ln x \right)^2 .$

Example: Consider the nonhomogeneous Euler equation of third order:

$x^3 y''' + x^2 y'' -2x\,y' +2\, y = x^3 \ln x \qquad (x>0) .$
The associated homogeneous equation $$x^3 y''' + x^2 y'' -2x\,y' +2\, y = 0$$ has three linearly independent solutions $$y_1 =x, \ y_2 = x^2 , \ y_3 = x^{-1} = 1/x$$ because the corresponding auxiliary polynomial $$P(m) = m\left( m-1 \right) \left( m-2 \right) + m \left( m-1 \right) -2m +2 = \left( m^2 -1 \right) \left( m-2 \right)$$ has three real nulls (where P(m) = 0): $$m_1 =1, \ m_2 = 2, \ m_3 = -1 .$$

According to variation of parameters, we seek a particular solution in the form:

$y_p (x) = A(x)\, x + B(x)\, x^2 + C(x) \, x^{-1} ,$
where unknown functions should be determined from the Lagrange system of equations for their derivatives:
\begin{align*} A' (x) \, x + B' (x) \,x^2 + C' (x) \, x^{-1} &= 0, \\ A' (x) + B' (x) \, 2x - C' (x) \, x^{-2} &= 0, \\ 2\, B' (x) + 2\, C' (x) \, x^{-3} &= \ln x . \end{align*}
We ask Mathematica to provide us a solution for the above system of algebraic equations, and it did not refuse:
Solve[{a*x + b*x^2 + c/x == 0, a + 2*b*x - c/x^2 == 0, 2*b + 2*c/x^3 == Log[x]}, {a, b, c}]
{{a -> -(1/2) x Log[x], b -> Log[x]/3, c -> 1/6 x^3 Log[x]}}
To find exact expressions of unknown functions A(x), B(x), and C(x), we just need to integrate:
\begin{align*} A' (x) &= - \frac{1}{2} \, x\, \ln x \qquad \Longrightarrow \qquad A(x) = \frac{x^2}{8} \left( 1 - 2\, \ln x \right) , \\ B' (x) &= \frac{1}{3}\, \ln x \qquad \Longrightarrow \qquad B(x) = \frac{1}{3}\left[ x\, \ln x -x \right] = \frac{x}{3}\left( \ln x -1 \right) , \\ C' (x) &= \frac{1}{6}\, x^3 \,\ln x \qquad \Longrightarrow \qquad C(x) = \frac{x^4}{96} \left( 4\,\ln x -1 \right) . \end{align*}
Plugging these expressions into the formula $$y_p (x) = A(x)\, x + B(x)\, x^2 + C(x) \, x^{-1} ,$$ we obtain a particular solution
$y_p (x) = A(x)\, x + B(x)\, x^2 + C(x) \, x^{-1} = \frac{x^3}{32} \left( 4\,\ln x - 7 \right) .$
■

1. Find a fundamental set of solutions of $$x^2 y'' + 2x\,y' + y = 0 .$$

1. Delkhosh, M., A method for solving the special type of Cauchy-Euler differential equations and its algorithms in MATLAB, 2012, Journal of Science, Vol. 2, No. 3, pp. 131--135.