function aa = mat1d(mtype, s, h, ix, x, b, aphib, xlen, blen, Lx)
%       Current version: June 27, 2017
%
%    We compute either the lower matrix (mtype = -1), the approximate
%    matrix (mtype = 0), or the upper matrix (mtype = 1) of size (xlen, xlen)
%
%    The computation is done by finding matrices corresponding to
%    values of b = m. The main issue here is to sum over the values
%    of b which give contributions from smallest to largest, to avoid 
%    loss of accuracy. 
%
%  We use the facts that for all x_1, x_2 in [0,1],
%  v_s(x_1) <= v_s(x_2) exp(2*s |x_1 -x_2|/gam),
%  2s(2s+1)/(2/gam  + Gam)^2 \le |v_s''(x)| \le 2s(2s+1)/gam^2 v_s(x)
%  where gam = b(blen), the smallest value in the vector b
%  and Gam = b(1), the largest value in the vector b.
%  From the error for linear interpolation, we have for
%  x_r <= theta_m(x) = 1/(x + m) <= x_{r+1},
%  [1 - err^1(x_k)] v_s^I(theta)m(x)) \le v_s(theta_m(x)) 
%  \le  [1 - err^2(x_k)]v_s^I(theta)m(x))
%  where err^1(x) = [x_{r+1} - theta_m(x)][theta_m(x) - x_r]
%  s (2s+1)/gam^2 exp(2sh/gam)
%  and err^2(x) = [x_{r+1} - theta_m(x)][theta_m(x) - x_r]
%  s (2s+1)/(2/gam + Gam)^2 exp(-2sh/gam)
%
     gam = b(blen);
     Gam = b(1);
%    Initialize matrix
aa = zeros(xlen, xlen);
for k = 1:xlen
  for l = 1:blen
    thetabxkl  = thetabx(x(k),b(l));
    ixkl = ix(k,l);
    v1 =  Lx(k,l,2);   v0 = Lx(k,l,1);
%   display([v0, v1])
%
%
%%%                   if (v0 < 0 | v0 >1)
%%%  display(['v0 < 0 | v0 >1', v0, thetabxkl, hint]) % error check
%%%                    end
%%%                    if (v1 < 0 | v1 >1) 
%%%  display(['v1 < 0 | v1 >1', v1, thetabxkl, hint]) % error check
%%%                    end
kk0 = ixkl;
kk1 = ixkl +1;
%  
%   The following statements check that we have not exceeded the array size
%
     if (kk0 > xlen) display(['kk0 > xlen', kk0, xlen, ixkl])
     end
     if (kk1 > xlen) display(['kk1 > xlen', kk1, xlen, ixkl])
     end
%
% 
if (mtype == -1)
   err = (s*(2*s+1)/gam^2)*exp(2*s*h/gam)*v0*v1*h*h;
elseif (mtype == 0) 
   err = 0.;
elseif (mtype ==1)
   err = (s*(2*s+1)/(2/gam + Gam)^2)*exp(-2*s*h/gam)*v0*v1*h*h;
else
  fprintf('mtype is not one of the allowed values');  
 end
%
    aphibpskl = aphib(k,l)^(2*s);
    aa(k,kk0) = aa(k,kk0) + aphibpskl*v0*(1.0-err);
    aa(k,kk1) = aa(k,kk1) + aphibpskl*v1*(1.0-err);
  end
end

% end of function mat1d


