Fast algorithm for the MATLAB's $\tt impz$ function? Does anybody here know what is the code algorithm behind the impz function? 
I know it uses the recursive algorithm but how is such implemented in impz that it can get the impulse response of very high $N$ and $M$ orders (numerator and denominator coefficients)? 


I can make a script that gets the impulse response of a second order LCCDE (linear constant coefficient differential equation) but I want to know a GENERAL code that can get the impulse response of any LCCDEs regardless of their order.
the impz function of matlab, simply computes N samples of the impulse response of a given LTI system from its LCCDE description coefficients a nd b.
As far as I know its efficiency is a secondary issue, since computing an impulse response in real time is not a frequently encountered task.

One equivalent code is simply this:

  % given a,b vectors and number of samples as N
  x = [1 zeros(1,N-1)];
  h = filter(b,a,x);


I guess you know how the filter function is implemented and what its output means. In the simplest terms, the filter() matlab function computes the output vector y for an input vector x of an LTI system: $ \sum_{k=0}^{k=N}a_k y[n-k] = \sum_{k=0}^{k=M} b_k x[n-k]$ , by the recursive approach: 

$$ y[n] = - \sum_{k=0}^{k=N} \frac{a_k}{a_0} y[n-k] + \sum_{k=0}^{k=M} \frac{b_k}{a_0} x[n-k] $$

Where the initial conditions of the system is assumed to be zero, inline with impulse response definition and LTI system with initial rest conditions. 

Note that in the impulse response computation, the signal $x[n]$ has only a single nonzero entry at $n=0$ and the rest of the computation is a recursion on the signal $y[n]$ alone.