Hay un cerrado fórmula para resolver este tipo de ecuaciones:
Deje $b,c,x_0,x_1$ ser números complejos, definir $x_2,x_3,\dots$ por
$$x_{n+2}+b\,x_{n+1}+c\,x_n=0,$$
deje $u$ $v$ ser la solución de $z^2+b\,z+c=0$, y poner
$$s_n:=u^n+u^{n-1}\,v+u^{n-2}\,v^2+\cdots+v^n$$
para $n\ge0$, que es
$$s_0:=1,\quad s_1:=u+v,\quad s_1:=u^2+uv+v^2,\dots$$
En particular
$$s_n=\frac{u^{n+1}-v^{n+1}}{u-v}\quad\text{if}\quad u\not=v$$
y
$$s_n=(n+1)\ u^n\quad\text{if}\quad u=v.$$
Entonces tenemos
$$x_n=s_{n-1}\ x_1-s_1\ s_{n-1}\ x_0+s_n\ x_0$$
para $n\ge1$.
Más generalmente, hay una (un poco menos explícita) cerrado fórmula para resolver ecuaciones de la forma
$$x_{n+k}+a_{k-1}\ x_{n+k-1}+\cdots+a_0\ x_n=b_n$$
en términos de $x_0,x_1,\dots,x_{k-1}$. Si alguien está interesado yo estaría encantado de darle.
EDIT. Gracias a Didier Piau, por su interés! [Añadido: ... y para señalar los errores de tal manera divertida!] La siguiente es una variación en torno a los casi obsceno y bien conocido el hecho de que una sola integración es suficiente para resolver el ODE $y^{(n)}=f$.
Deje $D$ ser un grado $q\ge1$ monic polinomio; deje $\mathbb C^\mathbb N$ el conjunto de secuencias de números complejos; deje $\Delta$ ser el operador de desplazamiento a la asignación de la secuencia de $x$ $\mathbb C^\mathbb N$ a la secuencia de $\Delta x$ $\mathbb C^\mathbb N$ definido por $(\Delta x)_t=x_{t+1}$; deje $f$$\mathbb C^\mathbb N$; deje $c_0,\dots,c_{q-1}$ ser números complejos; deje $y$ ser el único elemento de $\mathbb C^\mathbb N$ satisfactorio
$$D(\Delta)\,y=f,\quad y_n=c_n\text{ for all }n<q.$$
Vamos a dar
una fórmula para $y$ en términos del resto de la división larga de $X^t$$D$, y
una fórmula para el resto de cualquier tiempo la división de un polinomio $P$ por un polinomio distinto de cero $D$.
Para $(n,t)$ $\mathbb N^2$ denotar por $g_n(t)$ el coeficiente de $X^n$ en el resto de la división larga de $X^t$$D$. Aquí está la primera fórmula:
Si $t\ge q$ $\displaystyle y_t=\sum_{n<q}\,c_n\,g_n(t)+\sum_{k<t}\,g_{q-1}(t-1-k)\,f_k.$
La prueba de la primera fórmula. Introducir la secuencia de vectores $v_t:=(y_t,\dots,y_{t+q-1})$. Tenemos
$$v_{t+1}=B\,v_t+f_t\,e_q,\quad v_0=c,$$
donde $e_q$ es el último vector de la base canónica de $\mathbb C^q$, e $B$ es la transpuesta de la $q$ $q$ matriz $A$ caracteriza por
$$j<q\Rightarrow A\,e_j=e_{j+1},\quad A\,e_q=-a_0\,e_1-\dots-a_{q-1}\,e_q.$$
Por lo tanto
$$v_t=B^t c+f_0\,B^{t-1}\,e_q+f_1\,B^{t-2}\,e_q+\dots +f_{t-1}\,e_q.$$
Basta entonces para tomar el primer componente del lado izquierdo y derecho, y tener en cuenta el siguiente hecho:
Deje $D$ ser un grado $q\ge1$ polinomio, vamos a $P$ ser un polinomio, y
$$b_{q-1}\,X^{q-1}+\dots+b_0$$
el resto de la división larga de $P$$D$. Entonces tenemos
$$P(A)\,e_1=b_0\,e_1+\dots+b_{q-1}\,e_q.$$
QED
(2) Vamos a $R$ ser el resto de la división de un polinomio $P$ por un polinomio distinto de cero $D$. Aquí $P$$D$$\mathbb C[X]$. Podemos asumir
$$
D(X)=\big(X-a_1\big)^{m_1}\cdots\big(X-a_r\big)^{m_r},
$$
donde el $a_j$ son distintos y el $m_j$ positivo. Para cualquier fracción racional $f(X)\in\mathbb C(X)$ definido en $a_j$, vamos a $T_j(f(X))$ es el grado en que la mayoría de $m_j-1$ Taylor aproximación de $f(X)$$X=a_j$. La segunda fórmula es
$\displaystyle R(X)=\sum_{j=1}^r\,T_j\left(P(X)\,\frac{(X-a_j)^{m_j}}{D(X)}\right)
\frac{D(X)}{(X-a_j)^{m_j}}\quad.$
Si $m_j=1$ todos los $j$, obtenemos la Fórmula de Interpolación de Lagrange
$$
R(X)=\sum_{j=1}^r\,P(a_j)\,\prod_{k\=j}\,\frac{X-a_k}{a_j-a_k}\quad.
$$
La prueba de la segunda fórmula. Vamos a resolver las congruencias
$$R\equiv P\bmod(X-a_j)^{m_j}$$
como se explica en esta respuesta. QED