Here is the procedure.
Denote the two roots by $r_1$ and $r_2$, with $r_1 \gt r_2$.
The Method of Frobenius will always generate a solution corresponding to $r_1$, but may generate a solution for the smaller second root $r_2$ of the indicial equation.
If the method fails for $r_2$, then an approach is to keep the recursion solution in terms of $r$ and use it to find the coefficients $a_n$ (for $n \ge 1$), in terms of both $r$ and $a_0$, where $a_0 \ne 0$. For ease, $a_0$ is typically chose to be one.
Using this more general form and the coefficients, the two independent solutions can be written as:
$$y_1(r, x) = x^r \sum_{n=0}^\infty ~ a_n(r)x^n = \sum_{n=0}^\infty ~ a_n(r)x^{n+r} \\y_2(r,x) = \dfrac{\partial}{\partial r}[(r - r_2)y_1(r, x)]~\Bigr|_{r=r_2}$$
You should be able to use this approach and show:
$$y(x) = y_1(x) + y_2(x) = \dfrac{c_1(3x(3x^2 + 3x+ 2) + 2)}{x^3} + \dfrac{c_2e^{3x}}{x^3}$$
The differential equation
\begin{align}
y^{''} - x y^{'} + y = 0
\end{align}
can be solved via a power series of the form
\begin{align}
y(x) = \sum_{k=0}^{\infty} a_{n} x^{n} = a_{0} + a_{1} x + a_{2} x^{2} + \cdots .
\end{align}
It is fairly evident that
\begin{align}
\sum_{k=0}^{\infty} k(k-1) a_{k} x^{k} = \sum_{k=0}^{\infty} (k-1) a_{k} x^{k}
\end{align}
which yields the equation for the coefficients
\begin{align}
a_{k+2} = \frac{ (k-1) a_{k} }{ (k+1) (k+2) }.
\end{align}
It is discovered that $a_{3} = 0 \cdot a_{1}$. Since, for $k$ being odd, say $k \rightarrow 2k+1$,
\begin{align}
a_{2k+3} = \frac{k a_{2k+1} }{(k+1)(2k+3)}
\end{align}
it is clear that all the odd coefficients depend of $a_{3}$ for $k \geq 1$ and leads to $a_{2k+1} = 0$ for $k \geq 1$. The even $k$ values are
\begin{align}
a_{2} &= - \frac{a_{0}}{2!} \\
a_{4} &= - \frac{a_{0}}{4!} \\
a_{6} &= - \frac{(1 \cdot 3) a_{0}}{6!} \\
a_{8} &= - \frac{(1\cdot 3 \cdot 5)a_{0}}{8!}
\end{align}
which has the general form
\begin{align}
a_{2k} = - \frac{a_{0}}{2^{k} k! (2k-1)}.
\end{align}
The series for $y(x)$ now be seen in the form
\begin{align}
y(x) = a_{0} + a_{1} x - a_{0} \sum_{k=1}^{\infty} \frac{ x^{2n} }{2^{k} k! (2k-1)}.
\end{align}
The power series discovered can be evaluated as follows. Consider
\begin{align}
\partial_{x} \left( \sum_{k=1}^{\infty} \frac{ x^{2n-1} }{2^{k} k! (2k-1)} \right) &= \sum_{k=1}^{\infty} \frac{ x^{2n-2} }{2^{k} k!} = \frac{1}{x^{2}}( e^{x^{2}/2} -1).
\end{align}
Integrating both sides
\begin{align}
\sum_{k=1}^{\infty} \frac{ x^{2n} }{2^{k} k! (2k-1)} &= x \int^{x} \frac{e^{u^{2}/2} -1}{u^{2}} du = x \left[ \sqrt{\frac{\pi}{2} } erfi\left( \frac{x}{\sqrt{2}} \right) - \frac{e^{x^{2}/2}}{x} + \frac{1}{x} \right] \\
&= \sqrt{\frac{\pi}{2} } \cdot x \cdot erfi\left( \frac{x}{\sqrt{2}} \right) - e^{x^{2}/2} + 1.
\end{align}
With this series the general solution of $y(x)$ can be sen by
\begin{align}
y(x) &= a_{0} + a_{1} x - \sum_{k=1}^{\infty} \frac{ x^{2n} }{2^{k} k! (2k-1)} \\
&= a_{1} x - a_{0} \left[ \sqrt{\frac{\pi}{2} } \cdot x \cdot erfi\left( \frac{x}{\sqrt{2}} \right) - e^{x^{2}/2} \right],
\end{align}
where $erfi(x)$ is the imaginary error function (erfi(x) = -i erf(ix)).
Best Answer
This is not an answer but it is too long for a comment.
As commented by Winther, finding something nicer than the general recurence seem to be almost impossible.
Playing with a CAS and starting from what Winther wrote, what I obtained is the following monster for $a_{2n}$ $$\frac{ 4^{n-1} \cosh \left(\frac{\sqrt{7} \pi }{2}\right) \Gamma \left(\frac{1}{4} \left(1-i \sqrt{7}\right)\right) \Gamma \left(\frac{1}{4} \left(1+i \sqrt{7}\right)\right) \Gamma \left(n-\frac{i \sqrt{7}}{4}-\frac{1}{4}\right) \Gamma \left(n+\frac{i \sqrt{7}}{4}-\frac{1}{4}\right)}{\pi ^2 \Gamma (2 n+1)}a_0$$ $a_{2n+1}$ is as nice.
As Winther said, this kind of beauty just reflects the fact that the solution of the differential equation involves ugly hypergeomatric functions $$y=c_1 \, _2F_1\left(-\frac{1}{4}-\frac{i \sqrt{7}}{4},-\frac{1}{4}+\frac{i \sqrt{7}}{4};\frac{1}{2};x^2\right)+i c_2 x \, _2F_1\left(\frac{1}{4}-\frac{i \sqrt{7}}{4},\frac{1}{4}+\frac{i \sqrt{7}}{4};\frac{3}{2};x^2\right)$$