Hi Rik, your new mpoles() seems to work fine on my side. With your fix also residue() gives correct results in my case.

Many thanks!

This is actually two bugs in one.  The mpoles function needs to be fixed, but even with that in place residue also requires changes for the motivating example in this bug report.

I have modified mpoles.m, but while inspecting the code I found that it is very old (original check-in in 2007) and there are many improvements that should be made including using correct indentation, validating the inputs, having complete BIST tests, matching the variable names in the code to the documentation, and expanding the documentation to clearly state what the function does.  I've started on that.

Here is an example where residue returns a wrong result if at least one pole is zero:

p = [0 -1e8 -9e6 -8e5 -8e3];
r = [10 -1 -2 -3 -4];
[b, a] = residue(r, p, []);

In this case I get the wrong result
b =
   4.0000e+00   5.5963e+08   7.0093e+15   7.2560e+21   5.7597e+25

The correct result for b would be:
b =
            0   1.2043e+08   3.0605e+15   4.3760e+21   5.7600e+25

The result for 'a' is okay.

After a closer look I found that the problem is within the function 'mpoles':
If a pole is nonzero, the function regards it unique, if the others are separated by a factor larger than 'tol' (same as Matlab does). 
But if a pole is zero, it regards it unique, if the other poles have a magnitude larger than 'tol' times the mean of all poles.

The latter is actually the root cause for the problem. Matlab in contrast treats 'tol' as an absolute tolerance for poles equal zero.

Here is an example for 'mpoles':

p = [0 -1e8 -9e6 -8e5 -8e3];
tol = 1e-3; % Hard coded in 'residue'
[mults, i] = mpoles(p, tol, 0)

The output is:
mults = 1  2  1  1  1
i     = 1  5  2  3  4

The output from Matlab is:
mults = 1  1  1  1  1
i     = 1  2  3  4  5

Can you change the behavior of 'mpoles', so that is treats 'tol' as an absolute tolerance for poles equal zero, like Matlab does?