%STURAN_HERMITE Symbolic Gauss-Turan quadrasture formula
%  for Hermite weight function
%
%
global  s
format long e
hom=0;
s=1; N=4;
dig=64; eps0=.5e-62; eps1=.5e-14;
fprintf('digits=%3.0f,  s=%2.0f,  N=%3.0f\n',dig,s,N)
%
% select values of n for which complete reults 
% are to be printed
nprint=4;
%
R=vpa(zeros(N,2*N-1)); al0=vpa(zeros(5,1));
abn=vpa(zeros(N-1,2));
x0=[1.596071637983321;2.180607124035126;2.759002709962278; ...
3.333551485269060;3.905451757616834];
for j=1:5
  al0(j)=sroot(dig,@p,@pd,eps0,x0); 
end
R(1,1)=al0(s);
rho(1)=R(1,1);
sab0=sr_laguerre(dig,(s+1)*N);
fprintf('eps1=%12.4e\n',eps1)
for n=1:N
  print=nnz(nprint-n)<size(nprint,1);
  if n==1
    rho_init=rho(1);
  elseif n==2
    rho_init=[rho(1);vpa(linfun(2));vpa(quadfun(2))];
%
%   for the initial run (with N=4)
%
%    rho_init=[rho(1);3.7*sab0(2,1);11.3*sab0(2,2)];
%
  else
    rho_init(2*n-1)=vpa(quadfun(n));
%
%   for the initial run (with N=4)
%
%    if n==3, rho_init(2*n-1)=14.4*sab0(3,2); end
%    if n==4, rho_init(2*n-1)=15.3*sab0(4,2); end
%
    for k=n+1:2*n-2
      rho_init(k)=rho(k-1);
    end
    rho_init(n)=vpa(linfun(n));
%   for the initial run (with N=4)
%
%    if n==3, rho_init(n)=3.9*sab0(3,1); end
%    if n==4, rho_init(n)=3.9*sab0(4,1); end
%
    for k=1:n-1
      rho_init(k)=rho(k);
    end
  end
  fprintf('\n')
  tic
  [sxw,rho,it,condJmax,condAmax]=sturan(dig,n,s,eps1,sab0,rho_init,hom,print);
  time=toc;
  fprintf('n=%2.0f\n',n)
  fprintf('time=%8.2f\n',time)
  fprintf('it=%2.0f\n',it)
  fprintf('condJmax=%12.4e\n',condJmax)
  eps0=.5*10^(-(dig-log10(condJmax)-2));
  if eps0>eps1
    error('eps1 too small to cope with ill-conditioning')
  end
  R(n,1:2*n-1)=rho';
  if n>1
    abn(n-1,1)=R(n,n); abn(n-1,2)=R(n,2*n-1);
  end
%
% the constant terms in the linear resp. quadratic
% function determining initial values in Newton's method
%
%  if n==2 & N==4
%    fprintf('R(2,2)=%8.4f\n',subs(R(2,2)))
%    fprintf('R(2,3)=%8.4f\n',subs(R(2,3)))
%  end
%
  if print
    disp('                       alpha')   
    subs(R(n,1:n))'
    fprintf('\n')
    disp('                        beta')
    subs([0;R(n,n+1:2*n-1)'])
    fprintf('\n')
    for sig=1:2*s+2
      if sig==1
        disp('                       nodes')
      else
        disp('                      weights')
      end
      k=1:n;
      subs(sxw(:,sig))
      fprintf('\n')
    end
  end
end
%
% when N=4 determine linear and quadratic functions for
% computing initial approximations in Newton's method
%
%if N==4
%  a1=(R(4,4)-R(2,2))/2
%  b1=R(3,5)-R(2,3), d=R(4,7)-R(3,5); b2=(d-b1)/2
%
% the linear resp. quadratic functions are
%     R(2,2)+a1*(n-2),
%     R(2,3)+b1*(n-2)+b2*(n-2)*(n-3)
%
%end
n=(1:N-1)';
figure
hold on
plot(n,subs(abn(:,1)),'*')
n1=1; n2=N-1; 
a1=abn(n1,1); a2=abn(n2,1);
d=(a2-a1)/(n2-n1);
alin=a1+d*(n-n1);
plot(n,subs(alin),'-')
hold off
figure
hold on
plot(n,subs(abn(:,2)),'*')
n1=1; n2=ceil(N/2); n3=N-1;
a1=abn(n1,2); a2=abn(n2,2); a3=abn(n3,2);
d1=(a2-a1)/(n2-n1); d2=(a3-a2)/(n3-n2);
dd=(d2-d1)/(n3-n1);
bquad=a1+d1*(n-n1)+dd*(n-n1).*(n-n2);
plot(n,subs(bquad),'-')
hold off
