How to get the characteristic equation?

Here’s a rote rule for doing so. Start with the recurrence:

$$U_n=3U_{n-1}-U_{n-3}$$

Convert each subscript to an exponent:

$$U^n=3U^{n-1}-U^{n-3}$$

Change the variable to the one that you want to use in the characteristic equation:

$$x^n=3x^{n-1}-x^{n-3}$$

Divide through by the smallest exponent, in this case $n-3$:

$$x^{n-(n-3)}=3x^{(n-1)-(n-3)}-1\;,$$

which simplifies to $$x^3=3x^2-1\;.$$

With a little practice you can do the conversion in one go. For instance, the recurrence $$a_n=4a_{n-2}-6a_{n-3}+3a_{n-4}$$ has characteristic equation $$x^4=4x^2-6x+3\;,$$ as you can check by following through the steps given above.

Often one sees much handwaving here: rote rules, guesses, etc. But the conceptual background is very simple. Let $\rm\,S\,$ be the linear $\,n$-shift operator $\rm\:S\, u_n = u_{n+1}.\,$ In operator form we have

$$\begin{align} \rm 0\, &=\,\rm u_{\large n+3} - 3\, u_{\large n+2} + u_{\large n}\\[.2em] &=\,\rm S^{\large 3} u_{\large n} - 3\, S^{\large 2} u_{\large n} + u_{\large n}\\[.2em] &=\rm (\underbrace{S^{\large 3}\ \ -\ \ 3\,S^{\large 2} +\, 1}_{\rm\Large f(S)})\, u_{\large n}\\ \end{align}\quad$$

Factoring $\rm\,f(S) = (S-c_1)\,(S-c_2)\,(S-c_2)\,$ for $\rm\, c_j\in \mathbb C,\,$ reduces to solving linear recurrences $\rm\: (S-c)\,u_n = 0,\:$ i.e. $\rm\,u_{n+1} = c\,u_n,\,$ so $\rm\:u_n = u_o c^n.\,$ Because the $\rm\,c_j,\:$ are independent of $\rm\,n\,$ the operators $\rm\, S-c_j\:$ commute, so if the roots $\rm\,c_j\,$ are distinct, this yields three solutions $\rm\,c_j^n.\,$ Simple linear algebra (using the Casoratian) shows that these three solutions span the solution space.

Same for LDEs: $\rm\ D\!-\!ax\,$ kills $\,\rm e^{ax}$ so $\,\rm (D\!-\!i)(D\!+\!i) = D^2\!+\!1\,$ kills $\,\rm (e^{ix}\!+\!e^{-ix})/2=\cos(x)\,$

and $\rm\,S\!-\!x,\, S\!-\!y\,$ kill $\,\rm x^n,y^n$ so $\,\rm(S\!-\!x)(S\!-\!y) = S^2\! -\! (x\!+\!y) S\! +\! xy\,$ kills $\,\rm x^n+y^n =: f_n,\ $ i.e. $\rm \,f_{n+2}-(x\!+\!y)\,f_{n+1}+xy\, f_n = 0$.

Thus the characteristic polynomial is simply the polynomial $\rm\,f(S)\,$ or $\rm\,f(D)\,$ obtained from writing the difference / differential equation in operator form, and the form of the solutions follows immediately from factoring the characteristic polynomial. This reduces the solution of an $\rm\,n$-th order constant coefficient equation to the solution of $\rm\,n\,$ first-order constant coefficient equations.

Remark $\ $ The above depends crucially on the fact that the coefficients $\rm\,c_j\,$ are constant, i.e. independent of $\rm\,n,\,$ so they commute with $\rm\,S,\,$ i.e. $\rm\, Sc = cS\ $ (vs. $\rm\,Sc_{\large n} = c_{\large \color{#c00}{n+1}}S;\,$ $\rm Dc = c D+\color{#c00}{c'} 1)$. This property is what allows us to faithfully map the factorization of the characteristic polynomial $\rm\,f(x)\,$ into the same factorization for $\rm\,f(S),\, $ i.e. polynomial arithmetic assumes that $\,\rm x\,$ commutes with all coef's $\rm\,c_j,\,$ e.g $\rm\,x^2-c^2 = (x-c)(x+c)\iff cx = xc.\ $ The proof that the factorization of $\rm\,f(x)\,$ is true depends only on the ring axioms and the fact that $\rm\,x\,$ commutes with the coefficients, so the proof will persist in any ring where such constant commutativity persists. See this answer for further discussion from the viewpoint of the universal property of polynomial rings.

The same sort of sort of reduction works for any product of linear operators that commute.

"Guess" that $U(n) = x^n$ is a solution and plug into the recurrence relation:

$$ x^n = 3x^{n-1} - x^{n-3} $$

Divide both sides by $x^{n-3}$, assuming $x \ne 0$:

$$ x^3 = 3x^2 - 1 $$

Which is the characteristic equation you have.