iir filter saturation problems

I want to implement IIR filter without MATLAB functions. I get the coefficients as follows (Fs = 1000 , BPF 5-15Hz)

 f1=5;                                                                      
     f2=15;                                                                     
     Wn=[f1 f2]*2/fs;                                                           
     N = 3;                                                                     
     [b,a] = butter(N,Wn); 
     b = flipud(b')
     a = -flipud(a')

b coff:

1.0e-04 *
    
    -0.2915
          0
     0.8744
          0
    -0.8744
          0
     0.2915

a coff:

 -0.8819
         5.3942
       -13.7568
        18.7244
       -14.3455
         5.8657
        -1.0000

And I have written the filter code as below, but the output goes to infinity.

        for i = 7 : N
            BPF(i) =  0.0001*(-0.2915*X(i) +0.8744* X(i-2) - 0.8744* X(i-4) + 0.2915* X(i-6))...
                +(-0.8819*BPF(i-1)+5.3942*BPF(i-2)-13.7568* BPF(i-3)+18.7244*BPF(i-4)-14.3455*BPF(i-5)+5.8657*BPF(i-6));
            
        end

Answer 1

[Chunru]

It seems that your implementation of filter is not correct and hencie make the iir filter not stable:

for i = 7 : N
    BPF(i) =  0.0001*(-0.2915*X(i) +0.8744* X(i-2) - 0.8744* X(i-4) + 0.2915* X(i-6))...
        +(-0.8819*BPF(i-1)+5.3942*BPF(i-2)-13.7568* BPF(i-3)+18.7244*BPF(i-4)-14.3455*BPF(i-5)+5.8657*BPF(i-6));
    
end

Replace above with

BPF = filter(b, a, X) 

[nxp nxp]

its right.But why it is not stable? for initial value of output? I want to implement it with C

[Chunru]

For iir filter, the poles of transfer function must be within unit circle to be stable.

Your filter formula is implemting a different filter from "filter(b, a, x)" and it somehow not stable. Check out your implementation according the the IIR filter formula.