Difference Equation

M-order homogenous difference equation:
sigma(m = 0 : M - 1, a[m] * x[n+m]) = 0 (... 1.1)
can be represented by
sigma(m = 0 : M - 2, a⁽¹⁾[m] * x⁽¹⁾[n+m]) = 0 (... 1.2),
where x⁽¹⁾[k] = a"[1]*x[k+1] + a[0]*x[k]
a⁽¹⁾[k] = a"[k+1]/a"[1], k = 0 : M-1; a⁽¹⁾[M-2] = a[M-1]/a'[1]
a'[k], k = 1 : M-2, satisfy
a[M-1]/a'[M-2] = a"[1] / a[0] = a"[k+1]/a'[k], k = 1 : M-3
We suppose the solution for equation 1.2 is
x⁽¹⁾[n] = sigma(m = 0 : M - 2, b_mc_mⁿ)
say, a"[1]*x[n+1] + a[0]*x[n] = sigma(m = 0 : M - 2, b_mc_mⁿ)
and c_m is the solution set for characteristic equation: sigma(m = 0 : M - 2, a⁽¹⁾[m]*t^m) = 0 (... 1.3)
Let p = -a[0]/a"[1], we have x[n+1] = p * x[n] + sigma(m = 1 : M - 2, b_mc_mⁿ) (... 1.4).
The solution for the 2-order difference equation is:
x[n] = pⁿ * x[0] + sigma(m = 1 : M-2, b_m*sigma(k = 0 : n-1, c_m^kp_n^n-1-k)) (... 1.5)
After calculating the geometric series with ratio c_m/p, we can see the formula 1.5 is in the form
x[n] = sigma(m = 1 : M-2, b'_m*c_mⁿ)+b'_M-1*pⁿ
The characteristic equation for 1.1 is
sigma(m = 0 : M - 1, a[m] * t^m) = 0,
It can be rewritten as
sigma(m = 0 : M - 2, a⁽¹⁾[m] * t^m*(a"[1]*t+a[0])) = 0 (... 1.6)
So, c_m, m = 1 : M-2, plus p = -a[0]/a"[1] is the solution set for 1.6.
By this induction, the theorem that states the form of the solution for a homogenous difference equation is the weighted sum of the powers with degree equal to the index of the term, of the roots for the corresponding characteristic equation.

Another proof of this is in signal processing style:
First, we have to restrict the signal to be rational in this proof, say:
sigma(m=0:M-1, a[m]*x[n+m]) = 0, n>=0, and x[n] = 0, n < 0.
So, there exists sigma(n=0:M-1, a[m]*z^m*sigma(n=0:Inf, x[n+m]*z^-(n+m))) = 0.
sigma(n=0:Inf, x[n+m]*z^-(n+m)) = X(z) - sigma(n=0:m-1, x[n]*z^-n), where X(z) is the z-transform of sequence x.
From this we have X(z) = x[0] + Y(z) / sigma(m=0:M-1, a[m]*z^m) = x₀ + K(z)
The signal K(z) has poles: z_m, which are the roots of sigma(m=0:M-1, a[m]*z^m) = 0. Assured by the definition of the difference equation, it has no pole at origin, and by the fact that the order Y(z) is less than M-1,
we have K(z) in the form of partial fraction: K(z) = sigma(m=1:M-1, A_m / (z - z_m)).
1/(z-z_k) can be expanded as Laurent series around origin, within ROC z_m < |z| < Inf, which represents the rational part: X(z) = x[0] + sigma(n=1:Inf, sigma(m=1:M-1, A_mz_m^n-1)z^-n)
So the sequence is derived as: x[n] = sigma(m=1:M-1, A_mz_m^n-1), n>0. Proof done.

However, unfortunately, the theorem above is not friendly to computer algorithm design, bringingno benefit tofast and stable implementation, as far as I can see.
To calculate any specifiic stage of output of a linearsystem, another method on the basis of matrix is more feasible:

For anyM-order equation discussed above,it can berewritten as:
x[n+m] = sigma(k =0 : M-2, a[k] * x[n+k]), the matrix representative of whichis:
[O, Eye(M-2);
a[0 : M-2]]
e.g., for the extended Fibonacci sequence x[n+3] = x[n]+x[n+1]+x[n+2] (x[0]=x[1]=1,x[2]=2), the matrix is:
A = [ 0 1 0
0 0 1
1 11]
Any output with index greater than 2 is in the third row of
a = A^(n-2) * T, where T = [1,1,2]'
The temporal complexity of the algorithm is determined by that of the power of matrix, which is O(log(n)).
The algorithm outperforms the simple iterative ones only when n is large enough so that the overall advantagein complexity compensates the disadvantage in individual matrix multiplication.