%
% SR_FREUD  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 Freud polynomials orthogonal with respect to
%    the weight function
%           w(x)=|x|^alpha*exp(-|x|^beta).
%    The output variable ab is the Nx2 array of the 
%    nofdig-digit recurrence coefficients, and dig is 
%    the number of digits required to achieve this 
%    target precision
%
function [ab,dig]=sr_freud(N,alpha,beta,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=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);