% Heat and boundary equations in Cartesian coordinates 
function Domain = RecFDEqs(x,y,plate,qgen,q1,q2,q3,q4,varargin)
  if length(varargin)
    Domain = RecCellFDEqs(x,y,plate,qgen,q1,q2,q3,q4,varargin{1});
  else
    [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)alpha*((T(m,n-1)-2*T(m,n)+T(m,n+1))/dx^2 + ...
      (T(m-1,n)-2*T(m,n)+T(m+1,n))/dy^2+qgen(T,t,m,n)/k);         
    BN1=@(T,t,m,n)alpha*(2*(-q1(T,t,m,n)/k+(T(m,n+1)-T(m,n))/dx)/dx + ...
      (T(m-1,n)-2*T(m,n)+T(m+1,n))/dy^2+qgen(T,t,m,n)/k);         
    BN2=@(T,t,m,n)alpha*(2*(-q2(T,t,m,n)/k+(T(m,n-1)-T(m,n))/dx)/dx + ...
      (T(m-1,n)-2*T(m,n)+T(m+1,n))/dy^2+qgen(T,t,m,n)/k);         
    BN3=@(T,t,m,n)alpha*((T(m,n-1)-2*T(m,n)+T(m,n+1))/dx^2 + ...
     2*(-q3(T,t,m,n)/k+(T(m+1,n)-T(m,n))/dy)/dy+qgen(T,t,m,n)/k); 
    BN4=@(T,t,m,n)alpha*((T(m,n-1)-2*T(m,n)+T(m,n+1))/dx^2 + ...
     2*(-q4(T,t,m,n)/k+(T(m-1,n)-T(m,n))/dy)/dy+qgen(T,t,m,n)/k); 
    BN5=@(T,t,m,n)alpha*(2*(-q1(T,t,m,n)/k+(T(m,n+1)-T(m,n))/dx)/dx + ...
     2*(-q3(T,t,m,n)/k+(T(m+1,n)-T(m,n))/dy)/dy+qgen(T,t,m,n)/k); 
    BN6=@(T,t,m,n)alpha*(2*(-q2(T,t,m,n)/k+(T(m,n-1)-T(m,n))/dx)/dx + ...
     2*(-q3(T,t,m,n)/k+(T(m+1,n)-T(m,n))/dy)/dy+qgen(T,t,m,n)/k); 
    BN7=@(T,t,m,n)alpha*(2*(-q1(T,t,m,n)/k+(T(m,n+1)-T(m,n))/dx)/dx + ...
     2*(-q4(T,t,m,n)/k+(T(m-1,n)-T(m,n))/dy)/dy+qgen(T,t,m,n)/k); 
    BN8=@(T,t,m,n)alpha*(2*(-q2(T,t,m,n)/k+(T(m,n-1)-T(m,n))/dx)/dx + ...
     2*(-q4(T,t,m,n)/k+(T(m-1,n)-T(m,n))/dy)/dy+qgen(T,t,m,n)/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
end
