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, bringing no benefit to fast and stable implementation, as far as I can see. 
To calculate any specifiic stage of output of a linear system, another method on the basis of matrix is more feasible:

For any M-order equation discussed above, it can be rewritten as:
x[n+m] = sigma(k = 0 : M-2, a[k] * x[n+k]), the matrix representative of which is:
[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 1 1]
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 advantage in complexity compensates the disadvantage in individual matrix multiplication.