%STURAN_HERMITE_EXPLORE Symbolic Gauss-Turan quadrature formula
%  for Hermite weight function. Preliminary routine for purposes
%  of exploration.
%
global  s
format long e
hom=0;
s=1; N=16; 
dig=72; 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 results 
% are to be printed
nprint=0;
%
R=vpa(zeros(N,2*N-1)); 
x0=[1.053573782217666;1.609349223985799;2.165699416895354; ...
2.722281296487472;3.278979167347821];
be1=sroot(dig,@p,@pd,eps0,x0);
R(2,3)=be1;
rho=[R(2,1:2) be1]';
sab0=sr_hermite(dig,(s+1)*N);
fprintf('eps1=%12.4e\n',eps1)
for n=2:N
  print=nnz(nprint-n)<size(nprint,1);
  if n==2
    rho_init=rho;
  elseif n==3
    rho_init=[R(n,1:n) .9881*be1 vpa(linfunGTH(n))]';
%
%   for the initial run (with N=5)
%
%    rho_init=[R(n,1:n) .98561*be1 6.45945*sab0(n,2)]';
  else
    rho_init(2*n-1)=vpa(linfunGTH(n));
%
%   for the initial run (with N=5)
%
%    if n==4, rho_init(2*n-1)=6.39351*sab0(n,2); end
%    if n==5, rho_init(2*n-1)=6.34578*sab0(n,2); end
    for k=n+1:2*n-2
      rho_init(k)=rho(k-1);
    end
    rho_init(1:n)=vpa(0);
%
%   for the initial run (with N=5)
%
%    if n==4, rho_init(n)=vpa(0); end
%    if n==5, rho_init(n)=vpa(0); end
  end
%  subs(rho_init)
  fprintf('\n')
  tic
  [sxw,rho,it,Jmax,condJmax]=sturan(dig,n,s,eps1,sab0,rho_init,hom,print);
  time=toc;
  rho(1:n)=vpa(0);
%  subs(rho)
  fprintf('n=%2.0f  s=%2.0f\n',n,s)
  fprintf('time=%8.2f\n',time)
  fprintf('it=%2.0f\n',it)
  fprintf('Jmax=%12.4e\n',Jmax)
  fprintf('condJmax=%12.4e\n',condJmax)
  eps0=.5*10^(-(dig-log10(condJmax)-2));
  fprintf('eps0=%12.4e\n',eps0)
  if eps0>eps1
    error('eps1 too small to cope with ill-conditioning')
  end
  R(n,1:2*n-1)=rho';
  if print
    disp('                       alpha')   
    vpa(R(n,1:n)',24)
    fprintf('\n')
    disp('                        beta')
    vpa([0;R(n,n+1:2*n-1)'],24)
    fprintf('\n')
    for sig=1:2*s+2
      if sig==1
        disp('                       nodes')
      else
        disp('                      weights')
      end
      k=1:n;
      vpa(sxw(:,sig),24)
      fprintf('\n')
    end
  end
end
%error('forced break')
% when N=5 determine the linear function for computing
% initial approximations in Newton's method
%
%if N==5
%  db=(R(5,9)-R(2,3))/3
%
% the linear function is
%     R(2,3)+db*(n-2)
%
error('forced break')
C=zeros(N-1,2);
for k=1:N-1
  d=zeros(N-k,1);
  for m=k+1:N
    d(m-k)=subs(R(m,2*m-k));
  end
  if k<N-1
    d1=(d(N-k)-d(1))/(N-k-1); 
    C(k,1)=d1; C(k,2)=d(1)-d1;
    m=(1:N-k)';
    merrb(k)=max(abs((C(k,1)*m+C(k,2)-d)./d));
    if k==1
      figure
      plot(m,d,'*')
%      axis([0 N+1 0 N+1])
      hold on
      x=linspace(1,N)';
      plot(x,C(k,1)*x+C(k,2),'-')
      hold off
    end
  else
    C(N-1,1)=0; C(N-1,2)=R(N,N+1); 
  end
end
fprintf('n=%2.0f\n',n)
merrb'
subs(R)
C
