m=10;
theta=0.76;
for k=1:m
    for L=1:m
        Gaussiancopula=@(x,y) 1/(2.*pi.*(1-theta.^2).^(1/2)).*exp(-(x.^2-2.*theta.*x.*y+y.^2)/(2.*(1-theta.^2)));
        b=integral2(Gaussiancopula,-inf,norminv(k/m),-inf,norminv(L/m));
        if isnan(b)
           b=0;
        end
        c(k,L)=b;
    end
end