(PARI)
my(x='x, y='y, z='z, w='w);
R = 1/(1-(w*x*y*z+w*x*y+w*x*z+w*y+w*z+x*y+x*z+y+z));
diag(n, expr, var) = {
my(a = vector(n));
for (i = 1, #var, expr = taylor(expr, var[#var - i + 1], n));
for (k = 1, n, a[k] = expr;
for (i = 1, #var, a[k] = polcoeff(a[k], k-1)));
return(a);
};
diag(12, R, [x, y, z, w])