function [h,xim] = etaij(P,m,n,i,j)
h = 0;
xim = 0;
dims = repmat([m],[1 i]);
for xi = 1:m^i
  X = ind2coord(xi,dims);
  %hX = etaijX(P,m,n,i,j,X);
  hX = mxetaijX(P,m,n,i,j,X);  % might as well do it fast!
  %h = max(h,hX);
  if hX > h
    h = hX;
    xim = xi;
  end
end
return