function aa = mat1dpert(mtype, s, h, ix, x, b, aphib, xlen, blen, Lx, lam)

%  We compute the lower matrix (if mtype = -1) 
%  or the approximate matrix (if mtype =0) of size (xlen,xlen),
%  which is also an upper matrix.
%
%  We use the facts that for all x_1, x_2 in [0,1],
%  v_s(x_1) <= v_s(x_2) exp(K*s |x_1 -x_2|),
%  where K depends on lam and is computed below, and
%  |v_s''(x)| \le L v_s(x)
%  From the error for linear interpolation, we have for
%  x_r <= theta_m(x) <= x_{r+1},
%  [1 - err(x_k)] v_s^I(theta)m(x)) \le v_s(theta_m(x)) \le  v_s^I(theta_m(x))
%  where err(x) is computed below
%
%  Compute needed constants depending on mtype and lam
if (mtype == -1)
%   M0 = 35*lam/(4*(3+2*lam));
   if (lam <= 3/7) C1 = (2.5)*(3.5)*lam/(1 + 3.5*lam);
   else C1 = 3.5*lam*(3/(7*lam))^(3/5);
   end 
   if (lam <= 1/14) C2 = (15/4)*(7/2)*lam/(1 + 3.5*lam);
   else C2 = 3*(3.5*lam)^(4/5)/(4^(1/5));
   end
   E2 = 3.5*lam*5/(6 + 4*lam);
   C3 = (6+4*lam)/(4-3*lam);
   C4 = C1*C3;
   kap = (2 + 7*lam)/(6 + 4*lam);
%      err1 = s^2*C1^2*C3 + s*(6+4*lam)*C3*(C2 + C1*M0*C3)/(8+11*lam);
err1 = (s*C2 + 2*s^2*C1^2*kap/(1 - kap) + s*C1*E2/(1-kap))/(1 -kap^2);
end
%    Initialize matrices
      aa = zeros(xlen,xlen);
      for k = 1:xlen
       	 for l = 1:blen
         	   thetabxkl = thetabxlam(x(k),b(l),lam);
                   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;
           if (kk0 > xlen)
               display(['kk0 > xlen'])
               display([kk0, xlen, ixkl])  % error check
           end
           if (kk1 > xlen) 
               display(['kk1 > xlen', kk1, xlen, ixkl])  % error check
               display(['kk1 > xlen'])
               display([kk1, xlen, ixkl])  % error check
           end
%
if (mtype == -1)
      err = err1*exp(C4*s*h)*v0*v1*h*h;
elseif (mtype ==0)
   err = 0.0;
else
   fprintf('mtype is not one of the allowed values');  
end

           aphibskl = aphib(k,l)^(s);
           aa(k,kk0) = aa(k,kk0) + aphibskl*v0*(1.0 - err);
           aa(k,kk1) = aa(k,kk1) + aphibskl*v1*(1.0 - err);
         end
       end
%
% end of function mat1dpert
