% MOMBOSEEINSTEIN  Generates the first 2N moments to dig 
%    decimal places of the Bose-Einstein weight function 
%    w(x)=[x/(e^x-1)]^r on R_{+} with r a positive integer.
%
function mf=momboseeinstein(dig,N,r)
digits(dig); 
for k=1:2*N+r
  m(k,1)=gamma(vpa(k+1))*zeta(vpa(k+1));
  if r==1 & k<=2*N
    mf(k)=m(k,1);
  end
end
if r>1
  for rho=2:r
    for k=1:2*N+r-rho
      m(k,rho)=vpa(((k+rho-1)/(rho-1))*m(k,rho-1) ...
        -m(k+1,rho-1));
      if rho==r
        mf(k)=m(k,rho);
      end
    end
  end
end