function ab=sr_BE(i)
%SR_BE Recurrence coefficients for polynomials orthogonal with 
%respect to Bose-Einstein type weight functions.
%   SR_BE(I) generates the first 80 recurrence 
%   coefficients in 99-digit arithmetic for the
%   polynomials orthogonal with respect to the weight
%   function W(T)=T/(E^T-1) on [0,Inf] if I=1, in 
%   135-digit arithmetic for the polynomials orthogonal
%   with respect to the weight function W(T)=[T/(E^T-1)]^2
%   on [0,Inf} if I=2, and in 98-digit arithmetic for 
%   polynomials orthogonal with respect to the weight
%   function W(T)=T^{1/2}/(E^T-1) on[0,Inf] if I=3, and 
%   places them, rounded to 32 digits, into the 80x2 
%   array AB.
   
syms mom ab 
if i==1
  digits(99); dig=99;
  for k=1:160
    mom(k)=gamma(vpa(k+1,dig))*zeta(vpa(k+1,dig));
  end
elseif i==2
  digits(135); dig=135;
  for k=1:160
    mom(k)=gamma(vpa(k+2,dig))*(zeta(vpa(k+1,dig)) ...
    -zeta(vpa(k+2,dig)));
  end
elseif i==3
  digits(98); dig=98;
  for k=1:160
    mom(k)=gamma(vpa(k+1/2))*zeta(vpa(k+1/2));
  end
else
  error('i not in range')
end 
ab=schebyshev(dig,80,mom);
ab=vpa(ab,32);
