function [ab,dig]=sr_freud(N,alpha,beta,nofdig)
%SR_FREUD  Recurrence coefficients of the Freud orthogonal
%polynomials.
%   [AB,DIG]=SR_FREUD(N,ALPHA,BETA,NOFDIG) uses the classical
%   Chebyshev algorithm to compute the first N recurrence coefficients
%   to NOFDIG digits of the polynomials orthogonal with respect to
%   the weight function W(X)=|X|^ALPHA*EXP(-|X|^BETA) on [-INF,INF]
%   and places them into the Nx2 array AB. The output variable DIG
%   is the number of digits required to achieve the desired accuracy.

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=momfreud(dig,N,alpha,beta);
  if i==dig0
    ab0=schebyshev(dig,N,mom);
  else
    ab1=schebyshev(dig,N,mom);
    serr=vpa(abs(ab1-ab0),dig);
%    serr=vpa(abs((ab1-ab0)./ab1),dig);
%    serr=vpa(abs((ab1(:,2)-ab0(:,2))./ab1(:,2)),dig);
    err=subs(serr);
    maxerr=max(max(err)); 
    ab0=ab1;
   end
end
ab=vpa(ab1,nofdig);
