最后我自己实现了：

function eq = residue2sym( R,P,K )

syms s

eq = 0;

lastR = NaN;
lastP = NaN;
multiplicityCounter = 1;
for i = 1:length(R)
    rc = R(i);
    pc = P(i);

    if(~isnan(lastP) )
        if(pc == lastP && abs(lastR) < 1e-15)   % Quando existe multiplicidade, o R correspondente ao primeiro P
                                                % da multiplicidade deveria
                                                % ser 0, mas estranhamente
                                                % aparece como -4.4409e-16.
                                                % O que deveria ser uma
                                                % comparação com 0 então
                                                % vira uma comparação para
                                                % um módulo bem pequeno
            multiplicityCounter = multiplicityCounter + 1;
        else
            eq = eq + lastR/(s - lastP)^multiplicityCounter;
            multiplicityCounter = 1;
        end
    end
    lastR = rc;
    lastP = pc;
end
if(~isnan(lastP))
    eq = eq + lastR/(s - lastP)^multiplicityCounter;
end

 if(~isempty(K))
     eq = eq + sum(s.^(0:length(K)));
 end

eq = vpa( eq, 2);

end

字符串

function [F, r, p, k] = residue2sym(b, a)
    [r, p, k] = residue(b, a);
    syms s;
    repeat = ones(size(p));
    for i = 2:length(p)
        if p(i) == p(i-1)
            repeat(i) = repeat(i-1) + 1;
        end
    end
    F = sum(r./(s-p).^repeat) + poly2sym(k, s);
end

字符串