% SR_EXP0INFAB Symbolic/variable-precision recurrence
%    coefficients for the weight function 
%    w(x)=exp(-a/x-bx) on [0,inf]
%
syms a b sa sb c mom ab
N=80; a=1; b=1;
digits(120); dig=120; dstr=num2str(dig);
sa=vpa(a,dig); sb=vpa(b,dig);
c=vpa(2*sqrt(a*b),dig);
cstr=char(c);
maple(strcat('Digits:=',dstr));
K1str=strcat('evalf(BesselK(1,',cstr,'));');
K2str=strcat('evalf(BesselK(2,',cstr,'));');
m1str=maple(K1str);
m2str=maple(K2str);
m0=vpa(m1str,dig);
m1=vpa(m2str,dig);
mom(1)=vpa(2*sqrt(a/b)*m0,dig);
mom(2)=vpa(2*(a/b)*m1,dig);
for k=3:2*N
  mom(k)=((k-1)*mom(k-1)+sa*mom(k-2))/sb;
end
ab=schebyshev(dig,N,mom);
vpa(ab,40)
