function mua = eval2dPnse(Iset, s, N, J, ixp,ixm, iyp, iym, ...
aphibp, aphibm, zlent, blen,deg, Lxp, Lxm, Lyp, Lym, bn)

  %       Current version: October 4, 2016
%    We compute the approximate matrix and its largest eigenvalue mua.
%
%    For Iset = 3, 
%    the computation is done by finding matrices corresponding to
%    values of b = m + in first for which n >=0 (the matrix ap),
%    and then the matrix for which n <0, (the matrix am), and 
%    finally adding the two.  
%
% If Theta_b(z) \in R_{k,l} with lower left corner x_k,y_l, define the 
% interpolant v_s^I of v_s (of degree deg in each variable ) by:
% the usual tensor product Lagrange interpolant of a function 
% of two variables on a rectangle.
% We shall speak of the corners of the rectangles as the mesh points
% and the set of all points at which we interpolate as the
% interpolation points
% zlent is the dimension of the interpolation points
%    Initialize matrices
a = zeros(zlent, zlent);
ap = zeros(zlent, zlent);
am = zeros(zlent, zlent);

for k = 1:zlent
  for l = 1:blen
%
%   The following calculations use the function
%    kindexPn = kindxPn(il,jl,id,jd,deg,N,J)
%   residing in a matlab file of the same name. 
%  Suppose the interpolation points are ordered in each row from left to right
%  starting with the interpolation points at level y=0.
%  Given an interpolation point of the form 
%   ((il+(ideg-1)/deg)*h, (jl+ (jdeg-1)/deg)*h),
%  with il = N/2-J(jl+1), N/2 + J(jl+1), jl =  0, N/2, 
%  and 1 <= ideg <= deg+1, 1 <= jdeg <= deg+1,
%  we wish to determine the index k of that interpolation point 
%  in an ordered list beginning with 1 and ending with zlent
%
    aphibpskl = aphibp(k,l)^s;
    ixklp = ixp(k,l);     iyklp = iyp(k,l);
  for ideg = 1:deg+1
     for jdeg = 1:deg+1
	  kk(ideg,jdeg) = kindxPn(ixklp,iyklp, ideg-1, jdeg-1, deg,N,J);
%   The following statement checks that we have not exceeded the array size
%       if (kk(ideg,jdeg)  > zlent) 
%	 display([ideg, jdeg, kk(ideg,jdeg), ixklp, iyklp])
%       end
	   Lxyp(ideg,jdeg) = Lxp(k,l,ideg)*Lyp(k,l,jdeg);
ap(k,kk(ideg,jdeg)) = ap(k,kk(ideg,jdeg)) + aphibpskl*Lxyp(ideg,jdeg);
     end
  end

    if (bn(l) > 0)  & (Iset == 1 | Iset ==3)
%     display(bn(l))
    ixklm = ixm(k,l);     iyklm = iym(k,l);
    aphibmskl = aphibm(k,l)^s;
  for ideg = 1:deg+1
     for jdeg = 1:deg+1
	  kk(ideg,jdeg) = kindxPn(ixklm,iyklm, ideg-1, jdeg-1, deg,N,J);
%   The following statement checks that we have not exceeded the array size
%       if (kk(ideg,jdeg)  > zlent) 
%	 display('kk > zlent')
%        display  kk(ideg,jdeg)
%        display ([ideg,jdeg])
%         display ([ixklp, iyklp])
%         display (zlent)
%       end
	   Lxym(ideg,jdeg) = Lxm(k,l,ideg)*Lym(k,l,jdeg);
am(k,kk(ideg,jdeg)) = am(k,kk(ideg,jdeg)) + aphibmskl*Lxym(ideg,jdeg);
     end
   end
  end
 end
end
       a = ap + am;
%

% Computation of the largest eigenvalue of the matrix a
% by the Matlab routine eigs.  
mua = eigs(a,1);
% end of function eval2dPnse


