If you use the characteristic polynomial, it will indeed be a cubic, but it’s a cubic that factors very easily: it’s
$$x^3-3x_2+3x-1=(x-1)^3\;.$$
If you use generating functions, note that you can write the recurrence as
$$a_n=3a_{n-1}-3a_{n-2}+a_{n-3}+8-8[n=0]-7[n=1]-9[n=2]\;,\tag{1}$$
valid for all $n\in\Bbb Z$ if you make the blanket assumption that $a_n=0$ for all $n<0$. The last three terms are Iverson brackets and are there to make the recurrence yield the correct initial values.
Now multiply $(1)$ by $x^n$ and sum over $n\ge 0$:
$$\sum_{n\ge 0}a_nx^n=3\sum_{n\ge 0}a_{n-1}x^n-3\sum_{n\ge 0}a_{n-2}x^n+\sum_{n\ge 0}a_{n-3}x^n+8\sum_{n\ge 0}x^n-8-7x-9x^2\;.\tag{2}$$
If the generating function is $A(x)$, $(2)$ can be rewritten as
$$A(x)=3xA(x)-3x^2A(x)+x^3A(x)+\frac8{1-x}-8-7x-9x^2\;,$$
which you can then solve for $A(x)$ and resolve into partial fractions.
Linear Recurrence Equations have typical solutions $a_n=\lambda^n$. Using this, we can compute the possible values of $\lambda$ for this equation from
$$
\lambda^n=-\lambda^{n-1}+2\lambda^{n-2}
$$
which means, assuming $\lambda\ne0$, that
$$
\lambda^2+\lambda-2=0
$$
This is the characteristic polynomial for the recurrence
$$
a_n=-a_{n-2}+2a_{n-2}
$$
The characteristic polynomial is $x^2+x-2$ which has roots $1$ and $-2$. Thus, the sequence is $a_n=b(1)^n+c(-2)^n$. Plugging in the values for $n=0$ and $n=1$ gives
$$
a_n=3+8(-2)^n
$$
Best Answer
Let's prove by induction that:
$$a_n = \frac{3^n-1}{2}$$
For $n = 1$, we have $\frac{3^1-1}{2} = 1 = a_1$. So that's good so far.
Assume our formula holds for $1\leq k \leq n$. Then:
$$a_{n+1} = 3a_n + 1 = 3\frac{3^n-1}{2} + 1 = \frac{3^{n+1}-3+2}{2} = \frac{3^{n+1}-1}{2}$$
Wich is what we wanted to prove.
($a_n = 1 + 3 + 9 + \ldots + 3^{n-1}$. This is called geometric series. You can find quite a lot of info about it).