Determining the general solution of a third order ODE based on one known linearly independent solution I am trying to solve the following third-order ordinary differential equation with functional coefficients:
$$ {d^3f(x) \over dx^3} +
   \Bigg( e^{2x} - {4 \over 3} \Bigg) {df(x) \over dx} +
   {4 \over 3}\Bigg( e^{2x} + {4 \over 9} \Bigg) f(x) = 0 \tag1$$
I know that the general solution of the above equation consists of three linearly independent functions:
$$
f(x) = f_1(x) + f_2(x)+f_3(x) \tag2
$$
I also know that $f_1(x)=C_1e^{-{4\over 3}x}$ is one out of the three linearly independent functions.
My question is: Since I know $f_1(x)$, is there a way I can determine the other two linearly independent functions?
In my old college textbook, I found a theorem that states that if I know one particular solution of a second-order linear homogenous ordinary differential equation, I can find the other linearly independent solution by using the substitution:
$$ f_2(x)=g(x)f_1(x) \tag 3$$
where $g(x)$ satisfies the condition:
$$  ln\Bigg(\Bigg|{dg(x)\over dx}\Bigg|\Bigg) = 
    -2 ln\Bigg(\Bigg|f_1(x)\Bigg|\Bigg) - \int h(x) dx \tag4 $$
In the above equation, $h(x)$ is the functional coefficient in the general form of the linear homogeneous second-order differential equation:
$$ {d^2f(x) \over dx^2} + h(x) {df(x) \over dx}+ q(x)f(x) = 0 \tag 5$$
So I am hoping I can do something similar in the case of third-order differential equations.
 A: In most cases solvable linear ODE that are not scalar of order one are highly structures. Exhibiting the structure allows to reverse-engineer the construction of it and possibly get closer to solutions. Here one hint is the coefficient $e^{2x}$
$$
0=\left(D^3-\frac43D+\frac{16}{27}+e^{2x}(D+\frac43)\right)f(x)
=\left(D^2-\frac43D+\frac49+e^{2x}\right)\left(D+\frac43\right)f(x)
$$
Following this factorization, the equation can be split as
\begin{align}
f'(x)+\frac43f(x)&=g(x)\\
g''(x)-\frac43g'(x)+\left(\frac49+e^{2x}\right)g(x)&=0.
\end{align}
So if the reduced second equation can be solved for $g$, then its solution can be inserted into the solution formula for linear first-order DE to get $f$.
Or one can push the factor $e^{-\frac43x}$ from the left to $f$ to get
\begin{align}
0&=\left((D-\frac43)^2-\frac43(D-\frac43)+\frac49+e^{2x}\right)\left(D-\frac43+\frac43\right)e^{-\frac43x}f(x)
\\
&=\left(D^2-4D+4+e^{2x}\right)D\left(e^{-\frac43x}f(x)\right)
\end{align}
From the completed square in the first terms one can now further simplify by shifting $e^{2x}$ towards the middle
$$
0=\left(D^2+e^{2x}\right)e^{2x}D\left(e^{-\frac43x}f(x)\right)
$$
So if $h$ is a solution to $h''(x)+e^{2x}h(x)=0$, the solution $f$ can be reconstructed as $f(x)=e^{\frac43x}\int e^{-2x}h(x)\,dx$.
Parametrizing $h(x)=u(e^x)$ gives
$$
tu''(t)+u'(t)+tu(t)=0,
$$
which classifies somewhere in the realm of Bessel equations.
