function muaa = eval1dPn(s, h, ix, x, b, aphib, xlen, blen, deg)

%  We compute the approximate matrix aa(.,.) of size (xlen, xlen)
%  and its largest eigenvalue muaa.
%
%    Initialize matrices
      aa = zeros(xlen,xlen);
      for k = 1:xlen
       	 for l = 1:blen
             thetabxkl = thetabx(x(k),b(l));
             ixkl = ix(k,l);
             aphibskl = aphib(k,l)^(2*s);
             Delth = ones(deg+1,1);
             for i = 1:deg+1
		 kk(i) = ixkl + i-1;
		 if (kk(i) > xlen)
                    display(['kk(i) > xlen'])
	            display([i, kk(i), xlen, ixkl])  % error check
                 end
               	 for j = [1:i-1 i+1:deg+1]
	         Delth(i) = (thetabxkl - x(ixkl + j-1))/((i-j)*h/deg)*Delth(i);
                 end
		 aa(k,kk(i)) = aa(k,kk(i)) + aphibskl*Delth(i);
                 end
         end
       end
%
%  Computation of the largest eigenvalue of the matrix aa
%  by the Matlab routine eigs.
       muaa = eigs(aa,1);
%  end of function eval1dPn
