% EXAMPLE3_60 Computations for Example 3.60.
%
clear m mc mp iq idelta irout AB
global m mc mp iq idelta irout AB
f1='%10.0f %5.0f %16.4e\n';
f2='%10.0f %5.0f     negative nodes %2.0f\n';
f3='%10.0f %5.0f     complex nodes %3.0f\n';
fprintf('\n')
disp('         n     m         error       z')
mc=1; mp=0; iq=1; idelta=2; AB=[-Inf Inf];
N=20; eps0=10^4*eps; Mmax=100; irout=1;
for m=1:3
  [ab,Mcap,kount]=mcdis(N,eps0,@qherm,Mmax);
  for n=1:N
    jac=zeros(n,n);
    for i=1:n
      jac(i,i)=ab(i,1);
      if i>1, jac(i,i-1)=ab(i,2); end
      if i<n, jac(i,i+1)=1; end
    end 
    tau=sort(eig(jac));
    neg=find(tau<0&imag(tau)==0); zneg=nnz(neg);
    com=find(imag(tau)~=0); zcom=nnz(com)/2;
    nn=isempty(neg);
    nc=isempty(com);
    if nn&nc
      P=zeros(n);
      for nu=1:n
        P(1,nu)=1; 
        if n>1, P(2,nu)=tau(nu)-ab(1,1); end
        if n>2
          for mu=3:n
            P(mu,nu)=(tau(nu)-ab(mu-1,1))*P(mu-1,nu)...
              -ab(mu-1,2)*P(mu-2,nu);
          end
        end
      end
      e1=zeros(n,1); e1(1)=ab(1,2);
      lam=P\e1; ta=zeros(n,2);
      ta(:,1)=tau; ta(:,2)=lam./(tau.^(m+1));
      x=linspace(0,ta(n,1));
      errmax=max(abs(splinefunc(n,m,x,ta)-exp(-x.^2))); 
      fprintf(f1,n,m,errmax)
    elseif ~nn&nc
      fprintf(f2,n,m,zneg)
    elseif nn&~nc
      fprintf(f3,n,m,zcom)
    end
  end
  fprintf('\n')
end