% EXAMPLE3_52 Computations for Example 3.52.
%
fprintf('\n')
f1='%15.0f %12.3e %12.3e\n';
disp('              n       E2          Einf')
s=[2.4048255577 5.5200781103 8.6537279129];
N1=52; m=5; n=20;
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))'; p=zeros(100,n-m+1);
for k=1:5:n-m+1
  E2(k)=sqrt(sum(xw(:,2).*(qhat(:,k)-f0).^2));
  q=clenshaw(k-1,x,1,0,ab,c);
  p(:,k)=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;
  Einf(k)=max(abs(p(:,k)-bessel(0,x))); 
  fprintf(f1,k+m-1,E2(k),Einf(k))
end
