% FIGURE3_4 Computations for Figure 3.4.
%
hold on
s=[2.4048255577 5.5200781103 8.6537279129];
N1=52; m=5;
for n=5:7
  t=linspace(0,s(3),N1)';
  N=N1-2; 
  xw=zeros(N,2); d=ones(1,N);
  t0=t(2:N+1); f=bessel(0,t0);
  xw=[t0,((t0.^2).*(t0-s(1)).*(t0-s(2)).*(t0-s(3))).^2/N]; 
  ab=r_hahn(N-1,0,0);
  ab(:,1)=t0(1)+(t0(N)-t0(1))*ab(:,1)/(N-1); 
  ab(:,2)=((t0(N)-t0(1))/(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));
  ab1=ab; ab=chri7(N-4,ab1,0);
  ab1=ab; ab=chri7(N-5,ab1,0);
  f0=f./((t0.^2).*(t0-s(1)).*(t0-s(2)).*(t0-s(3)))+(1+...
    (1/s(1)+1/s(2)+1/s(3))*t0)./(s(1)*s(2)*s(3)*(t0.^2));
  [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=1-(x./s(1)).^2+((s(1)+s(2))/(s(1)^2*s(2)^2))*...
  (x.^2).*(x-s(1))-(s(1)*s(2)+s(1)*s(3)+s(2)*s(3))*...
  (x.^2).*(x-s(1)).*(x-s(2))/(s(1)^2*s(2)^2*s(3)^2)+...
  (x.^2).*(x-s(1)).*(x-s(2)).*(x-s(3)).*q;
  bess=bessel(0,x);
  if n==5
    plot(x,bess) 
    plot(x,p,'-.')
  elseif n==6
    plot(x,p,'--')
  else
    plot(x,p,'.')
  end
end  
axis([0 9 -.5 1.1])
plot([0 9],[0 0])
xlabel('x')
ylabel('Bessel')
hold off
