Error when using residue function

Please tell me if there's something with my code or it just because of the function "residue"

          syms z
          dz = (2*z + 1)^4*(3*z + 1)^2*(z + 1)^3;
      % get all coefficients of the above polynomial, dz (largest degree --> smallest one)
          cs = fliplr( double(coeffs(fz)) );
      % partial fraction decomposition of "1/dz" using "residue" function
          [r,p,k]  = residue(1,cs);

The result shows that

    p = [-1.0000 -1.0000 -1.0000 -0.5005 -0.5005 -0.4995 -0.4995 -0.3333 -0.3333]

However, should it be like this to be correct? (notice the differences of the 4th-7th elements in two results)

    p = [-1.0000 -1.0000 -1.0000 -0.5000 -0.5000 -0.5000 -0.5000 -0.3333 -0.3333]

Also by using r,p,k achieved from "residue" function, I cannot convert the partial fraction expansion back to polynomial coefficients. How can I fix it?

[Jan]

Please explain, why you assume that the second value of p is correct.

[Luan Nguyen Thanh]

I think there are 3 reasons:

1) As explained in "https://www.mathworks.com/help/matlab/ref/residue.html#outputarg_p", p represents the poles of 1/dz, so I think it should be

p = [-1.0000 -1.0000 -1.0000 -0.5000 -0.5000 -0.5000 -0.5000 -0.3333 -0.3333]

2) Same as 1), you can see that -0.5005 and -0.4995 can not be the poles of 1/dz.

3) I explained in my question that with the coefficients "r,p,k" obtained from using residue, I cannot recover the partial fraction expansion to polynomial coefficients, which means when I write:

[r,p,k] = residue(1,cs)
[b2,a2] = residue(r,p,k)

I will get b2 and a2 different from 1 and cs, respectively.

Answer 1

[Star Strider]

If you have R2015a or later, use the partfrac (link) function.