Is there a widely available implementation of an adaptive, recursive, numerically stable IIR filter? I need to identify the coefficients of a linear, causal, time-invariant physical system that can be described by a classical state-space formulation.

For the sake of the example, suppose that the system has two states, only one of which is observed, and one input. Suppose that this state-space is described by the following matrices:

  a = 
               x1          x2
   x1   -1.25e-05    1.25e-05
   x2   1.389e-05  -2.083e-05

  b = 
              u1
   x1          0
   x2  6.944e-06

  c = 
       x1  x2
   y1   1   0

  d = 
       u1
   y1   0


The discrete, equivalent, zeroth-order-hold transfer function is given by:

  3.481e-05 z + 3.446e-05
  -----------------------
   z^2 - 1.97 z + 0.9704


In filter terminology, this is clearly an IIR filter because of the higher-order coefficients in the denominator.

I'm looking for a suitable adaptive filter algorithm that will:


be numerically stable
allow for online updates, i.e. be recursive


From my research so far, I've read that if you treat the recorded signal values as inputs to this filter, then it's possible to treat this filter as a FIR filter. This is called the equation error approach. One can then apply, for example, the QR-RLS adaptive FIR filter to find the coefficients.

But in my trials on the real system I often find filter coefficients that are not physically possible, and I just don't know if my hand-crafted adaptive filter is correct or not. I would like to try with an "official" implementation of such an adaptive filter.

Therefore, my question is whether there exists (for MATLAB or in C) a freely available implementation of an adaptive IIR filter that has the same advantages of a QR-RLS FIR adaptive filter.
While I may not be the final authority, I have looked into this. In short, no. 

From my own research, the problem is stability. When the filter is extended to have a wider bandwidth, which of course means that the y[n] values and lags 'move quickly', and then is quenched to a lower bandwidth, where the y[n] values and lags 'move slowly', the y[n] buffer may be loaded with fast oscillated entries that appear as initial conditions to the low bandwidth filter. The filter can easily go unstable. If the bandwidth is altered in an adiabatic manner this instability may be avoided (I'm supposing) but then there's an additional constraint on your adaptivity. 

Now, you don't have to take my word for it. In "Adaptive Filters: Theory and Applications" (2nd) page 5 (Amazon), Boroujeny writes, "However, as we shall see in the later chapters, because of the many difficulties involved in adaptation of IIR filters, their application in the area of adaptive filters is rather limited." The author states that the adaptation process can place (digital) poles outside of the unit circle, even if they started within the unit circle. The rest of his book is mostly dedicated to FIR filters.

I hope this information helps.