function outvar = PBRsysFullModel(v,y,FeedN2,FeedH2,FeedAr,FeedCH4 ...
    ,k0f,Ef,k0r,Er,eps,nu,G0,Dp,Ar,MWN2,MWH2,MWNH3,MWAr,MWCH4,CpCoef)
global Rgas Rgas2;
% unpack dependent variables
MflowN2 = y(1);
MflowH2 = y(2);
MflowNH3 = y(3);
P = y(4);
T = y(5);
H = y(6);
% total molar flow
TotFlow = MflowN2 + MflowH2 + MflowNH3 + FeedAr + FeedCH4;
% partial pressures in atm
PN2 = MflowN2/TotFlow*P;
PH2 = MflowH2/TotFlow*P;
PNH3 = MflowNH3/TotFlow*P;
% forward and reverse reaction rates
rf = k0f*exp(-Ef/Rgas/(T+273.15))*PN2^0.5*PH2^1.5;
rr = k0r*exp(-Er/Rgas/(T+273.15))*PNH3;
outvar(1) = -(rf-rr)*eps;  % differential balance on N2 in kgmol/h/m3
% algebraic equations from stoichiometry for H2 and NH3
outvar(2) = MflowH2 - ( FeedH2 - 3*( FeedN2 - MflowN2) );
outvar(3) = MflowNH3 - 2*( FeedN2 - MflowN2 );
% average molecular weight
MWavg = (MflowN2*MWN2 + MflowH2*MWH2 + MflowNH3*MWNH3 + FeedAr*MWAr ...
    + FeedCH4*MWCH4) / TotFlow;
% gas density
z = [ MflowN2 ; MflowH2 ; MflowNH3 ; FeedAr ; FeedCH4 ]/TotFlow;
w = [ 0.039 ; -0.216 ; 0.25 ; 0.001 ; 0.011 ];
Tc = [ 126.2 ; 33.19 ; 405.65 ; 150.86 ; 190.564 ];  % K
Pc = [ 3394 ; 1297 ; 11277 ; 4870 ; 4641 ];  % kPa
K = [ 0       -0.036  0.222  0       0.036 ;
      -0.036  0       0      0       0.202 ;
      0.222   0       0      0       0 ;
      0       0       0      0       0.023 ;
      0.036   0.202   0      0.023   0 ];
SV = SpecVol(T,P,z,w,Tc,Pc,K) ;
RhoGas = MWavg / SV ; % kg/m3
% gas viscosity from kinematic viscosity and density
mu = nu * RhoGas; % Pa-s
% differential Ergun equation, dP/dV in atm/m3
outvar(4) = -(150*(1-eps)/Dp/G0*mu+7/4)*(1-eps)/eps^3/Dp/RhoGas*G0^2 ...
    /Ar/101325;
HtRx = HtRxnP(T,P);  % heat of reaction
outvar(5) = T - findT(H,T,MflowN2,MflowH2,MflowNH3,FeedAr,FeedCH4,CpCoef);
% differential energy balance
outvar(6) = (rf-rr)*(-HtRx)*eps;
outvar = outvar';


