I have obtained the following z-domain transfer function:
$\frac{Y(z)}{U(z)}=\frac{3.3641×10^{-7}×z^{6}+1.5584×10^{-5}×z^{5}+6.7263×10^{-5}×z^{4}+ 5.5016×10^{-5}×z^{3}+8.525×10^{-6}×z^{2}+1.2303×10^{-7}×z+9.7492×10^{-20}}{z^{7}-5.5998×z^{6} +13.5548×z^{5}-18.3685×z^{4}+15.0393×z^3-7.4352×z^2+2.0541×z-0.2245}$
This is a 7-th order Butterworth filter with a corner frequency of 498.79 Hz.
Now, to obtain the difference equation so I can implement it on a microcontroller, I am multiplying the right side by $\frac{z^{-7}}{z^{-7}}.$
I thus obtain:
$\frac{Y(z)}{U(z)}=\frac{3.3641×10^{-7}×z^{-1}+1.5584×10^{-5}×z^{-2}+6.7263×10^{-5}×z^{-3}+
5.5016×10^{-5}×z^{-4}+8.525×10^{-6}×z^{-5}+1.2303×10^{-7}×z^{-6}+9.7492×10^{-20}×z^{-7}}{1-5.5998×z^{-1}
+13.5548×z^{-2}-18.3685×z^{-3}+15.0393×z^{-4}-7.4352×z^{-5}+2.0541×z^{-6}-0.2245×z^{-7}}$
I have done this because it allows me to work with values from the past, rather than try and read the future (both input and output). After re-arranging, I get the difference equation:
$y[n]=5.5598×y[n-1]-13.5548×y[n-2]+18.3685×y[n-3]-15.0393×y[n-4]+7.4352×y[n-5]-2.0541×y[n-6]+0.2245×y[n-7]+3.3641×10^{-7}×u[n-1]+1.5584 ×10^{-5}×u[n-2]+6.7263×10^{-5}×u[n-3]+5.5016×10^{-5}×u[n-4]+8.525×10^{-6}×u[n-5]+1.2303×10^{-7}×u[n-6]+9.7492×10^{-20}×u[n-7]$
I implemented this difference equation, but I am not getting the same results as for its z-domain transfer function. I couldn't identify any mistake in what I did. Could someone please help me identify it? I am 90% sure there is some mistake in my workings somewhere. I tried to do the same procedure for another (simpler) transfer function and the behaviour of the discrete difference equation matched exactly that of the simulated z-domain transfer function for the same input.
UPDATE 1:
The poles of the z-domain are: $0.8894±0.2805j, 0.7979±0.1995j, 0.747±0.1022j, 0.731$ and they are all stable.
Here's more proof it's stable (z domain):
I first constructed the following continuous transfer function, which I used together with the MATLAB c2d() function to get the z-domain transfer function I mentioned earliler. The method was "impulse" and a sampling frequency of 10 kHz. The continuous form is:
$\frac{Y(s)}{U(s)}=\frac{2.9696×10^{24}}{s^{7}+1.4084×10^{4}×s^6 +9.9181×10^7×s^5+4.4917×10^{11}×s^4+1.4077×10^{15}×s^3+3.053×10^{18}×s^2+ 4.2582×10^{21}×s+2.9696×10^{24}}$