Matlab for Method 2

... as in Method 1 for poles ...

% Construct a resonator as a sum of arbitrary modes 
% with unit residues, adding a near-zero at DC.

B = zeros(1,2*nsec+1);
impulse = [1,zeros(1,2*nsec)];
for i=1:nsec, 
    % filter() does polynomial division here:
    B = B + filter(A,denoms(i,:),impulse);
end;

% add a near-zero at DC 
B = filter([1 -0.995],1,B);

... as in Method 1 for display ...