% Heat and boundary equations in cylindrical coordinates 
function Domain = CylFDEqs(r,z,plate,qgen,q1,q2,q3,q4,varargin)
  if length(varargin)
    Domain = CylCellFDEqs(r,z,plate,qgen,q1,q2,q3,q4,varargin{1});
  else
    [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)alpha*((T(m,n-1)-2*T(m,n)+T(m,n+1))/dr^2+...
      (T(m,n+1)-T(m,n-1))./rr(m,n)/2/dr+(T(m-1,n)-...
      2*T(m,n)+T(m+1,n))/dz^2+qgen(T,t,m,n)/k);                     
    BN1=@(T,t,m,n)alpha*(2*(-r1*q1(T,t,m,n)/k+(r1+dr/2)*(T(m,n+1)-...
      T(m,n))/dr)/dr/(r1+dr/4)+(T(m-1,n)-...
      2*T(m,n)+T(m+1,n))/dz^2+qgen(T,t,m,n)/k);                     
    BN2=@(T,t,m,n)alpha*(2*(-rN*q2(T,t,m,n)/k+...
      (rN-dr/2)*(T(m,n-1)-T(m,n))/dr)/dr/(rN-dr/4)+(T(m-1,n)-...
      2*T(m,n)+T(m+1,n))/dz^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))/dr^2+...
      (T(m,n+1)-T(m,n-1))./rr(m,n)/2/dr+2*(-q3(T,t,m,n)/k+...
      (T(m+1,n)-T(m,n))/dz)/dz+qgen(T,t,m,n)/k);                    
    BN4=@(T,t,m,n)alpha*((T(m,n-1)-2*T(m,n)+T(m,n+1))/dr^2+...
      (T(m,n+1)-T(m,n-1))./rr(m,n)/2/dr+2*(-q4(T,t,m,n)/k+...
      (T(m-1,n)-T(m,n))/dz)/dz+qgen(T,t,m,n)/k);                    
    BN5=@(T,t,m,n)alpha*(2*(-r1*q1(T,t,m,n)/k+(r1+dr/2)*(T(m,n+1)-...
      T(m,n))/dr)/dr/(r1+dr/4)+2*(-q3(T,t,m,n)/k+(T(m+1,n)-...
      T(m,n))/dz)/dz+qgen(T,t,m,n)/k);                              
    BN6=@(T,t,m,n)alpha*(2*(-rN*q2(T,t,m,n)/k+(rN-dr/2)*(T(m,n-1)-...
      T(m,n))/dr)/dr/(rN-dr/4)+2*(-q3(T,t,m,n)/k+(T(m+1,n)-...
      T(m,n))/dz)/dz+qgen(T,t,m,n)/k);                              
    BN7=@(T,t,m,n)alpha*(2*(-r1*q1(T,t,m,n)/k+(r1+dr/2)*(T(m,n+1)-...
      T(m,n))/dr)/dr/(r1+dr/4)+2*(-q4(T,t,m,n)/k+(T(m-1,n)-...
      T(m,n))/dz)/dz+qgen(T,t,m,n)/k);                              
    BN8=@(T,t,m,n)alpha*(2*(-rN*q2(T,t,m,n)/k+(rN-dr/2)*(T(m,n-1)-...
      T(m,n))/dr)/dr/(rN-dr/4)+2*(-q4(T,t,m,n)/k+(T(m-1,n)-...
      T(m,n))/dz)/dz+qgen(T,t,m,n)/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
end
