Numerical Errors in mldivide

I wrote the following code to perform 2D heat transfer FEA with convection from the surface and quadrilateral elements. Currently it is a 2x2 element grid and heat is being generated in the first element. I would expect the result to be symmetric given the conduction in x and y is the same. However, the results are not symmetric. I believe I am assembling the K and F matrices correctly. Could the asymmetry be due to numerical errors with mldivide?

main file:

 clc; clear all; close all;

 eX = 2; %# X elements
 eY = 2; %# Y elements
 nodesPerElement = 4; %Quadrilateral element

 kx = 10; %x conductivity
 ky = 10; %y conductivity
 L = 0.05; %length of element
 W = 0.05; %width of element
 h = 1; %convective heat transfer coefficient
 z = 0.01; %thickness of board
 Tinf = 25; %ambient temperature

 graph = meshGrid(eX,eY);

numE = eX*eY; %# of elements
numN = (eX+1)*(eY+1); %# of nodes

Kglobal = zeros(numN,numN);
Fglobal = zeros(numN,1);

for e = 1:numE
    if e == 1
        qv = 1000;
    else
        qv = 0;
    end

      Kcx = kx*L/(6*W)*[-2 2 1 -1;
          2 -2 -1 1;
          1 -1 -2 2
          -1 1 2 -2];

      Kcy = ky*W/(6*L)*[-2 -1 1 2;
          -1 -2 2 1
          1 2 -2 -1
          2 1 -1 -2];

      Kh = -L*W*h/(36*z)*[4 2 1 2;
          2 4 2 1;
          1 2 4 2;
          2 1 2 4];

      K = Kcx + Kcy + Kh;
      F = -(h/z*Tinf + qv)*L*W/4*[1; 1; 1; 1];

      for n = 1:nodesPerElement
          for m = 1:nodesPerElement
              Kglobal( graph(e,n), graph(e,m)) = Kglobal( graph(e,n), graph(e,m)) + K(n,m);
          end

          Fglobal(graph(e,n)) = Fglobal(graph(e,n)) + F(1);
      end
  end

T = full(sparse(Kglobal)\sparse(Fglobal));

T = transpose(reshape(T, eX+1, eY+1));
x = linspace(0,eX*W,eX+1);
y = linspace(0,eY*L,eY+1);

[X,Y] = meshgrid(x,y);
surf(X,Y,T);
xlabel('x');
ylabel('y');
zlabel('Temperature (^oC)')
colorbar
axis equal

meshGrid.m:

function [map] = msehGrid(eX, eY)

numElements = eX*eY;
map = zeros(numElements, 4);

for n = 1:eY
    for m = 1:eX
        e = (n-1)*eX+m;
        map(e,1) = m + (n-1)*(eX+1);
        map(e,2) = m + 1 + (n-1)*(eX+1);
        map(e,3) = m  + (n)*(eX+1);
        map(e,4) = m + (n)*(eX+1) + 1;
    end
end

[Walter Roberson]

"Could the asymmetry be due to numerical errors with mldivide?"

Yes, if you count standard floating point round-off limitations as "numerical errors".