For situations involving linear, homogenous recurrence relations, the characteristic polynomial method works best.
Let us have the recurrence relation, for constants $c_i$ and $k>0$,
$$a_n = c_1 a_{n-1} + c_2 a_{n-2} + ... + c_{n} a_{n-k}$$
Much like with ordinary differential equations, we develop a characteristic polynomial for $a_n$. Let $a_n = \alpha^n$. Then we get
$$\alpha^n = c_1 \alpha^{n-1} + c_2 \alpha^{n-2} + ... + c_{k} \alpha^{n-k}$$
Divide through by $\alpha^{n-k}$ next:
$$\alpha^{k} = c_1 \alpha^{k-1} + c_2 \alpha^{k-2} + ... + c_{k}$$
Bring everything to the same side:
$$\alpha^{k} - c_1 \alpha^{k-1} - c_2 \alpha^{k-2} - ... - c_{k} = 0$$
This is a polynomial, and we seek the roots to it. Let them be $\alpha_1,...,\alpha_k$. (If $\alpha_i$ is a duplicate root, replace the first duplicate with $n\alpha_i$, the second with $n^2 \alpha_k$, and so on.)
Then, for these roots, up to constants $A_1,...,A_k$ depending on your initial conditions, we have
$$a_n = A_1 \alpha_1^n + ... + A_k \alpha_k^n$$
Footnotes & Caveats:
As you might imagine, it is hypothetically possible for the recurrence to have complex roots. If you're anything like me, this might present some worry especially for examples grounded in reality. You shouldn't worry too much about these cases. In a sense it all sort of "works out." Since you know one complex number is a root of the polynomial, so is its conjugate. You'll be raising those conjugates to the same power, and thus they remain conjugate. That is to say, for a complex number $z$,
$$\overline{z^n} = (\overline z)^n$$
Thus the conjugation remains and a lot of the "nasty" (imaginary) bits just so happen to negate as a result! A similar thing happens in the example I show in a bit, albeit not so much "complex conjugates" as just the conjugate of a rooted expression.
(A huge thanks for antkam in the comments for clearing up this detail and resolving a long-standing lack of an intuitive understanding as to how it all "works out" for me!)
Also, a nice tidbit: it's a good paranoia check to double-check your solution. Once you have the explicit form for $a_n$, check that your initial solutions are valid, and perhaps a few other values you obtain from the recurrence relation. In examples like this where you have to derive the recurrence relation yourself instead of simply being given it, you should be able to get some values by brute force for $n=1,2,3,$ and so on, for ever-how-many initial conditions you need to use. (You need as many initial conditions as there are previous values that determine $a_n$.)
Also bear in mind that this method only works for linear, homogenous recurrence relations. For nonhomogenous ones, I've spoken to you on solving them. For nonlinear ones, we need something more elaborate (such as generating functions) but such discussion is well beyond the scope of this post.
A simple example to motivate this method:
Example: Let us find the solution to the Fibonacci recurrence
$$a_n = a_{n-1} + a_{n-2}$$
where $a_0 = 0,$ and $a_1 = 1$.
(Bear in mind that while here each $a_{\text{something}}$ has coefficient $1$, they need not be, and as in the previous explanation the coefficients "carry over" to the characteristic polynomial. The Fibonacci relation is simply a common first example.)
Here, the characteristic polynomial is given by
$$\alpha^n = \alpha^{n-1} + \alpha^{n-2}$$
Divide through by $\alpha^{n-2}$:
$$\alpha^2 = \alpha + 1 \implies \alpha^2 - \alpha - 1 = 0$$
This quadratic can be solved by the quadratic formula. It's a well known result that the two roots to this are the golden ratio $\varphi$ and its conjugate $\overline \varphi$:
$$\varphi = \frac{1 + \sqrt 5}{2} \;\;\;\;\; \overline \varphi = \frac{1 - \sqrt 5}{2}$$
Thus, up to constants $A_1,A_2$, we can claim
$$a_n = A_1 \varphi ^n + A_2 \overline{\varphi}^n = A_1 \left( \frac{1 + \sqrt 5}{2} \right)^n + A_2 \left( \frac{1 - \sqrt 5}{2} \right)^n$$
What remains is to determine the constants $A_1, A_2$. To do this, substitute your initial conditions. Thus, you get a system of equations as below. In $a_0$, you let $n=0$ in your solution for $a_n$ above; similarly, $n=1$ in the $a_1$ case.
$$\left\{\begin{matrix}
a_0 = 0\\
a_1 = 1
\end{matrix}\right. \implies \left\{\begin{matrix}
A_1 + A_2 = 0\\
A_1 \varphi + A_2 \overline{\varphi} = 1
\end{matrix}\right.$$
To solve this is a fairly typical exercise in solving systems of equations, or linear algebra if you're faced with the awful situation of many initial conditions. I'll skip the boring bits, leaving the algebra to you, simply saying you should get $A_1 = 1/\sqrt 5, A_2 = -1/\sqrt5$.
And thus we get a general formula for the Fibonacci relation!
$$a_n = \frac{ \varphi ^n}{\sqrt 5} - \frac{\overline{\varphi}^n}{\sqrt 5} = \frac{ 1}{\sqrt 5}\left(\frac{1 + \sqrt 5}{2}\right)^n - \frac{1}{\sqrt 5}\left( \frac{1 - \sqrt 5}{2} \right)^n$$
