%MOD_MD  Modification by multiple divisors
%
function ab=mod_md(N,m,z,ab0,rho0)
for mu=1:m
  r0(mu)=feval(rho0,z(mu));
end
r(1)=r0(1); 
if m>1
  s=1;
  for mu=2:m
    s=-s;
    k=1:mu;
    for nu=1:mu
      c(nu)=1/prod(z(find(k-nu))-z(nu));
    end
    r(mu)=s*sum(c.*r0(1:mu));
  end 
end
ab=rschri4(N,m,ab0,z,r);
