Root locus for an open loop transfer function?
Ignoring the semantics of "root locus," you can certainly plot the roots of a polynomial as its coefficients change.  The following script plots the set of poles of your system as \$RL\$ varies from 0 to 10000 with \$IL\$ fixed at 1:

IL = 1;
RLvec = linspace(0,10e3,1e5);
rts = zeros(1e5,4);
for k = 1:length(RLvec);
  RL = RLvec(k);
  D = [8*IL, 28*IL+8*RL+3, 4*(7*IL+7*RL+3), 4*(2*IL+7*RL+3), 8*RL+4];
  r = roots(D);
  rts(k,:) = r';
end;
figure
hold on;
for b = 1:4;
  plot(real(rts(:,b)),imag(rts(:,b)));
end;
xlim([-5,0]); grid on;


It executes for me in under 5 seconds.

Admittedly, it does not show how the roots move as both parameters change, but that task is complicated by having a two-dimensional domain and a two-(real)-dimensional range.