%  Heat and boundary (cell) equations in cylindrical coordinates 
function Domain = CylCellFDEqs(r,z,plate,qgen,q1,q2,q3,q4,j)
  [alpha,k]=plate{1:2};
  dr=abs(diff(r(1:2))); dz=abs(diff(z(1:2))); % find dr and dz
  r1=r(1); rN=r(end);                         % inner and outer radius of domain
  rr=repmat(r,length(z),1);                   % creates 2D array rr(m,n)
  HEq=@(T,t,m,n,j)alpha*((T{j}(m,n-1)-2*T{j}(m,n)+T{j}(m,n+1))/dr^2+...
    (T{j}(m,n+1)-T{j}(m,n-1))./rr(m,n)/2/dr+(T{j}(m-1,n)-...
    2*T{j}(m,n)+T{j}(m+1,n))/dz^2+qgen(T,t,m,n,j)/k);
  BN1=@(T,t,m,n,j)alpha*(2*(-r1*q1(T,t,m,n,j)/k+(r1+dr/2)*(T{j}(m,n+1)-...
    T{j}(m,n))/dr)/dr/(r1+dr/4)+(T{j}(m-1,n)-...
    2*T{j}(m,n)+T{j}(m+1,n))/dz^2+qgen(T,t,m,n,j)/k);
  BN2=@(T,t,m,n,j)alpha*(2*(-rN*q2(T,t,m,n,j)/k+...
    (rN-dr/2)*(T{j}(m,n-1)-T{j}(m,n))/dr)/dr/(rN-dr/4)+(T{j}(m-1,n)-...
    2*T{j}(m,n)+T{j}(m+1,n))/dz^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))/dr^2+...
    (T{j}(m,n+1)-T{j}(m,n-1))./rr(m,n)/2/dr+2*(-q3(T,t,m,n,j)/k+...
    (T{j}(m+1,n)-T{j}(m,n))/dz)/dz+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))/dr^2+...
    (T{j}(m,n+1)-T{j}(m,n-1))./rr(m,n)/2/dr+2*(-q4(T,t,m,n,j)/k+...
    (T{j}(m-1,n)-T{j}(m,n))/dz)/dz+qgen(T,t,m,n,j)/k);
  BN5=@(T,t,m,n,j)alpha*(2*(-r1*q1(T,t,m,n,j)/k+(r1+dr/2)*(T{j}(m,n+1)-...
    T{j}(m,n))/dr)/dr/(r1+dr/4)+2*(-q3(T,t,m,n,j)/k+(T{j}(m+1,n)-...
    T{j}(m,n))/dz)/dz+qgen(T,t,m,n,j)/k);
  BN6=@(T,t,m,n,j)alpha*(2*(-rN*q2(T,t,m,n,j)/k+(rN-dr/2)*(T{j}(m,n-1)-...
    T{j}(m,n))/dr)/dr/(rN-dr/4)+2*(-q3(T,t,m,n,j)/k+(T{j}(m+1,n)-...
    T{j}(m,n))/dz)/dz+qgen(T,t,m,n,j)/k);
  BN7=@(T,t,m,n,j)alpha*(2*(-r1*q1(T,t,m,n,j)/k+(r1+dr/2)*(T{j}(m,n+1)-...
    T{j}(m,n))/dr)/dr/(r1+dr/4)+2*(-q4(T,t,m,n,j)/k+(T{j}(m-1,n)-...
    T{j}(m,n))/dz)/dz+qgen(T,t,m,n,j)/k);
  BN8=@(T,t,m,n,j)alpha*(2*(-rN*q2(T,t,m,n,j)/k+(rN-dr/2)*(T{j}(m,n-1)-...
    T{j}(m,n))/dr)/dr/(rN-dr/4)+2*(-q4(T,t,m,n,j)/k+(T{j}(m-1,n)-...
    T{j}(m,n))/dz)/dz+qgen(T,t,m,n,j)/k);
  FDEqs={BN1;BN2;BN3;BN4;BN5;BN6;BN7;BN8;HEq};
  M=length(z); N=length(r);
  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,r,z,k/alpha};
end
