%  Heat and boundary (cell) equations in Cartesian coordinates 
function Domain = RecCellFDEqs(x,y,plate,qgen,q1,q2,q3,q4,j)
  [alpha,k]=plate{1:2};
  if length(plate) == 3
    w=plate{end};
  else w=1;
  end
  dx=abs(diff(x(1:2))); dy=abs(diff(y(1:2))); % find dx and dy
  HEq=@(T,t,m,n,j)alpha*((T{j}(m,n-1)-2*T{j}(m,n)+T{j}(m,n+1))/dx^2 + ...
    (T{j}(m-1,n)-2*T{j}(m,n)+T{j}(m+1,n))/dy^2+qgen(T,t,m,n,j)/k);
  BN1=@(T,t,m,n,j)alpha*(2*(-q1(T,t,m,n,j)/k+(T{j}(m,n+1)-T{j}(m,n))/dx)/dx + ...
    (T{j}(m-1,n)-2*T{j}(m,n)+T{j}(m+1,n))/dy^2+qgen(T,t,m,n,j)/k);
  BN2=@(T,t,m,n,j)alpha*(2*(-q2(T,t,m,n,j)/k+(T{j}(m,n-1)-T{j}(m,n))/dx)/dx + ...
    (T{j}(m-1,n)-2*T{j}(m,n)+T{j}(m+1,n))/dy^2+qgen(T,t,m,n,j)/k);
  BN3=@(T,t,m,n,j)alpha*((T{j}(m,n-1)-2*T{j}(m,n)+T{j}(m,n+1))/dx^2 + ...
   2*(-q3(T,t,m,n,j)/k+(T{j}(m+1,n)-T{j}(m,n))/dy)/dy+qgen(T,t,m,n,j)/k);
  BN4=@(T,t,m,n,j)alpha*((T{j}(m,n-1)-2*T{j}(m,n)+T{j}(m,n+1))/dx^2 + ...
   2*(-q4(T,t,m,n,j)/k+(T{j}(m-1,n)-T{j}(m,n))/dy)/dy+qgen(T,t,m,n,j)/k);
  BN5=@(T,t,m,n,j)alpha*(2*(-q1(T,t,m,n,j)/k+(T{j}(m,n+1)-T{j}(m,n))/dx)/dx + ...
   2*(-q3(T,t,m,n,j)/k+(T{j}(m+1,n)-T{j}(m,n))/dy)/dy+qgen(T,t,m,n,j)/k);
  BN6=@(T,t,m,n,j)alpha*(2*(-q2(T,t,m,n,j)/k+(T{j}(m,n-1)-T{j}(m,n))/dx)/dx + ...
   2*(-q3(T,t,m,n,j)/k+(T{j}(m+1,n)-T{j}(m,n))/dy)/dy+qgen(T,t,m,n,j)/k);
  BN7=@(T,t,m,n,j)alpha*(2*(-q1(T,t,m,n,j)/k+(T{j}(m,n+1)-T{j}(m,n))/dx)/dx + ...
   2*(-q4(T,t,m,n,j)/k+(T{j}(m-1,n)-T{j}(m,n))/dy)/dy+qgen(T,t,m,n,j)/k);
  BN8=@(T,t,m,n,j)alpha*(2*(-q2(T,t,m,n,j)/k+(T{j}(m,n-1)-T{j}(m,n))/dx)/dx + ...
   2*(-q4(T,t,m,n,j)/k+(T{j}(m-1,n)-T{j}(m,n))/dy)/dy+qgen(T,t,m,n,j)/k);
  FDEqs={BN1;BN2;BN3;BN4;BN5;BN6;BN7;BN8;HEq};
  M=length(y); N=length(x);
  Nodes={2:M-1,1;2:M-1,N;1,2:N-1;M,2:N-1;1,1;1,N;M,1;M,N;2:M-1,2:N-1};
  Domain={FDEqs,Nodes,x,y,k/alpha*w};
end
