Solving two simultaneous recurrence relations

If we have the two recurrence relations $$a_n = 3a_{n-1} + 2b_{n-1}$$ $$b_n = a_{n-1} + 2b_{n-1}$$ with $a_0 = 1$ and $b_0 = 2$.

My solution is that we first add two equations and assume that $f_n = a_n + b_n$. The result is $f_n = 4f_{n-1}$. This can be solved easily and the solution is $f_n = a_n + b_n = 4^nf_0=4^n(3)$. We subtract the two equations from each other and we get $a_n - b_n = 2a_{n-1}$. By adding these two results we get $2a_n = 2a_{n-1} +4^n (3) $.

$a_n = a_{n-1} + (3/2) 4^n$ is annihilated by $(E-1)(E-4)$. Thus, the Generic solution is $a_n = \sigma + \beta4^n$. The constants $\sigma, \beta$ satisfy the equations: $a_0 = 1 = \sigma + \beta$ and $a_1 = 7 = \sigma + 4\beta$. Thus, the final solution is $a_n = -1 + (2) 4^n$ and $b_n = 4^n + 1$. Is this right? Thank you!

The procedure is right. There is no real need to wonder about correctness, since substituting the answers arrived at shows that they do satisfy the system.

Alternately, we can express the recurrence in matrix form, and use tools from linear algebra to find the solution.

For a lower level manipulation that works well, we can subtract as you did to get $$a_n-b_n=2a_{n-1},\quad\text{or equivalently}\quad b_n=a_n-2a_{n-1}.\tag{$1$}$$

Bump the indices in the first given recurrence to get $a_{n+1}=3a_n+2b_n$. Now use $(1)$ to substitute for $b_n$, and we get a second-order recurrence for the sequence $(a_n)$.

Remark: The lower level manipulation that we did does not depend on happy numerical "accidents." If we have a system $a_n=pa_{n-1}+qb_{n-1}$, $b_n=ra_{n-1}+sb_{n-1}$, use the two equations to eliminate say $b_{n-1}$, bump the indices as we did, and substitute.

Use generating functions, $A(z) = \sum_{n \ge 0} a_n z^n$ and similarly $B(z)$. Write: $$ \begin{align*} a_{n + 1} &= 3 a_n + 2 b_n \\ b_{n + 1} &= a_n + 2 b_n \end{align*} $$ By properties of generating functions: $$ \begin{align*} \frac{A(z) - a_0}{z} &= 3 A(z) + 2 B(z) \\ \frac{B(z) - b_0}{z} &= A(z) + 2 B(z) \end{align*} $$ With the given initial values: $$ \begin{align*} A(z) &= \frac{1 + 3 z}{1 - 4 z + z^2} \\ B(z) &= \frac{2 - 5 z}{1 - 4 z + z^2} \\ \end{align*} $$ Next step is to split $A(z)$ and $B(z)$ into partial fractions, and expand the results by geometric series to extract the terms. Sadly, in this case the zeros of the denominator are surds, and the expŕessions are unwidely.