How to solve two dimensional integral equation using gaussian quadratures in MATLAB

I need help with writing matlab code to solve two dimensional integral equation using gaussian quadratures. The following code solves one dimensional integral equation using gaussian quadratures % u(x) = 1 + 1/2.* int_{0}^{pi/4} sec^{2} (x) u(t) dt. n=5; [t,w]=lgwt(n,0,pi/4);% generate the quadrature points and weights for i=1:n for j=1:n k(i,j)=w(j).*0.5.*(sec(t(i))).^2; omk(i,j)=kronecker(i,j)-k(i,j); end g(i)=1; end f=omk\(g');

I tried to extend it to two dimensional case but did not get the correct answer.

for example to solve %% f(x,t) - int_{-1}^{1} int_{-1}^{1} (x sin(t) + y exp(s)) f(y,s) dy ds = x exp(-t) + 4x sin(t) - 7/3. I wrote the code n=32; a1 = -1; b1 = 1; a2 = -1; b2 = 1; [c1,w1]=lgwt(n,a1,b1); [c2,w2]=lgwt(n,a2,b2); A=eye(n); for l = 1:n for m=1:n

        for i = 1:n
            for j = 1:n
                K(i,j) = ( c1(l).*sin(c2(m)) + c1(i).*exp(c2(j)) ).*w1 (i).*w2 (j) ;
                A(i,j)=A(i,j) - K(i,j) ;
            end
        end
        G(l,m)= c1(l).*exp(-c2(m)) + 4.*c1(l).*sin(c2(m)) - 7/3;
    end
end
f = A\G'
But did not get the exact result.
 Any help will be appreciated.