% FIGURE3_3 Computations for Figure 3.3.
%
hold on
s=[2.4048255577 5.5200781103 8.6537279129];
N1=51; m=3; 
for n=3:5
  t=linspace(0,s(3),N1)';
  N=N1-1; 
  xw=zeros(N,2); d=ones(1,N);
  t0=t(1:N); f=bessel(0,t0);
  xw=[t0,((t0-s(1)).*(t0-s(2)).*(t0-s(3))).^2/N]; 
  ab=r_hahn(N-1,0,0);
  ab(:,1)=t0(N)*ab(:,1)/(N-1); ab(:,2)=(t0(N)/(N-1))^2*ab(:,2);
  ab(1,2)=1;
  ab1=ab; ab=chri7(N-1,ab1,s(1));
  ab1=ab; ab=chri7(N-2,ab1,s(2));
  ab1=ab; ab=chri7(N-3,ab1,s(3));
  f0=f./((t0-s(1)).*(t0-s(2)).*(t0-s(3)));
  [qhat,c]=least_squares(n-m,f0,xw,ab,d);
  x=(linspace(0,s(3)))'; 
  q=clenshaw(n-m,x,1,0,ab,c);
  p=(x-s(1)).*(x-s(2)).*(x-s(3)).*q;
  bess=bessel(0,x);
  if n==3
    plot(x,bess) 
    plot(x,p,'-.')
  elseif n==4
    plot(x,p,'--')
  else
    plot(x,p,'.')
  end
end  
plot([0 9],[0 0])
xlabel('x')
ylabel('Bessel')
hold off