% SR_FERMIDIRAC  This computes the first N recurrence 
%    coefficients alpha_k, beta_k, k=0,1,\ldots,N-1, to an 
%    accuracy of nofdig digits for the system of Fermi- 
%    Dirac polynomials orthogonal with respect to the 
%    weight function 
%                 w(x)=[1/(e^x+1)]^r. 
%    The oputput variable ab is the Nx2 array of the nofdig- 
%    digit recurrence coefficients, and dig is the number 
%    of digits required to achieve the target precision of 
%    nofdig decimal places.
%
function [ab,dig]=sr_fermidirac(N,r,nofdig)
syms mom ab ab0 ab1 
dd=10; dig0=nofdig;
i=dig0-dd; 
maxerr=1;
while maxerr>.5*10^(-nofdig)
  i=i+dd; dig=i;
  mom=momfermidirac(dig,N,r);
  if i==dig0
    ab0=schebyshev(dig,N,mom);
  else
    ab1=schebyshev(dig,N,mom);
    serr=vpa(abs((ab1-ab0)./ab1),dig);
    err=subs(serr);
    maxerr=max(max(err)); 
    ab0=ab1;
  end
end
ab=vpa(ab1,nofdig);