Your approach will work fine for the difference equation. $$k(n+1)3^{n+1}=k(n)3^{n+1}+2^{n-1}-1$$
Hence $3^{n+1}(k(n+1)-k(n))=2^{n-1}-1$, which you can solve the same way as in the differential equation case. In discrete-land $f(n+1)-f(n)$ serves the same role that $\frac{d}{dx}f(x)$ does in continuous-land.
Nb. In neither case would I choose this approach to use; it's usually easiest to use undetermined coefficients with a suitable guess, as @Qiaochu suggests.
Note: I changed the terminology somewhat; this sequence starts with $a_0$ rather than $a_1$.
Suppose we have a sequence $a_0,a_1,a_2,\dots$ whose generating function is
$$
f_a(x)=a_0+a_1x+a_2x^2+\cdots
$$
satisfying the recurrence relation
$$
a_n-P\,a_{n-1}-Q\,a_{n-2}-R\,a_{n-3}=0\iff\\
a_n=P\,a_{n-1}+Q\,a_{n-2}+R\,a_{n-3}
$$
Multiply $f_a(x)$ by the polynomial $1-Px-Qx^2-Rx^3$ to get the polynomial
$$
g(x) = b_0+b_1 x+ b_2 x^2+\cdots
$$
where for $n\geq 3$, $b_n=a_n-P\,a_{n-1}-Q\,a_{n-2}-R\,a_{n-3}$. By our recurrence relation, this means that $b_n=0$ whenever $n\geq 3$. So, we have
$$
(1-Px-Qx^2-Rx^3)f_a(x)=b_0+b_1 x+ b_2 x^2
$$
Which is to say that
$$
f_a(x)=\frac{b_0+b_1 x+ b_2 x^2}{1-Px-Qx^2-Rx^3}
$$
Where
$$
b_0 = a_0\\
b_1 = a_1 - P\,a_0\\
b_2 = a_2 - P\,a_1 - Q\,a_0
$$
Can you take it from there?
So in order to bring this back to the characteristic equation, we just need to use another little trick. Instead of writing this as a function of $x$, write it as a function of $\frac1x$. You could do this by making a substitution like $x=\frac1\omega$, but I prefer a more direct approach.
We have:
$$
f_a(x)=\frac{b_0+b_1 x+ b_2 x^2}{1-Px-Qx^2-Rx^3}
$$
With $b_1,b_2,b_3$ as defined above. From there, just divide the top and bottom by $x^3$ to get
$$
f_a(x)=\frac{b_0\left(\frac1{x}\right)^3+b_1 \left(\frac1{x}\right)^2
+ b_2 \left(\frac1{x}\right)}{
\left(\frac1{x}\right)^3-P\left(\frac1{x}\right)^2-
Q\left(\frac1{x}\right)-R}
$$
Now, suppose we have one repeated root. That is, $t^3 - Pt^2 - Qt - R=(t-r_1)(t-r_2)^2$ for roots $r_1,r_2$. We then can write the above as
$$
f_a(x)=\frac{b_0\left(\frac1{x}\right)^3+b_1 \left(\frac1{x}\right)^2
+ b_2 \left(\frac1{x}\right)}{
\left(\left(\frac1{x}\right)-r_1\right)
\left(\left(\frac1{x}\right)-r_2\right)^2}
$$
Where would you go from there? For the case of a triply repeated root, we have $t^3 - Pt^2 - Qt - R=(t-r)^3$ for the repeated root $r$. We then can write the generating function as
$$
f_a(x)=\frac{b_0\left(\frac1{x}\right)^3+b_1 \left(\frac1{x}\right)^2
+ b_2 \left(\frac1{x}\right)}{
\left(\left(\frac1{x}\right)-r\right)^3}
$$
Where would you go from there?
Best Answer
There is no general method (that i am aware of) , unless the coefficients in your linear recurrence relations are constant (which is not the case in your example).
Concerning your specific example, since $2n=(n+2)+(n-2)$, it rewrites $(n+2)b_{n+1}=(n-2)b_n$, where $b_n=a_{n+1}-a_n$. You can easily find $b_n$ from the former relation, then find $a_n$ since $a_n-a_0=\displaystyle\sum_{k=0}^{n-1}b_k$.