function r=dRr_dy(y,T,P)
% y : Mole fractions      [mol/mol]
% P : Reactor pressure    [bar]
% T : Reactor temperature [K]

global NC

I=ones(size(y))';
r=zeros(4,1);
H2=1;
Ga=2;
Pa=3;
Na=4;
Ar=5;

P=P*100;               % Total pressure [kPa]
p=P*y;                 % Component Partial pressure [kpa]
k=kf(T);
K=Keq(T);
dpdy=P*eye(NC);  

r(1,1:NC)=k(1)/K(1)*(K(1)*dpdy(Na,:)-dpdy(Ar,:)*p(H2)^3-3*p(Ar)*p(H2)^2*dpdy(H2,:)); % N <=> A + 3H2  [mol/(kg*s)]   Conversion of naphthenes to aromatics
r(2,1:NC)=k(2)/K(2)*(K(2)*(dpdy(Na,:)*p(H2)+p(Na)*dpdy(H2,:)) -dpdy(Pa,:));          % N + H2 <=> P   [mol/(kg*s)]   Conversion of naphthenes to paraffins
r(3,1:NC)=k(3)*dpdy(Na,:)/P;                                                         % N + H2 <=> 2G  [mol/(kg*s)]   Hydrocracking of naphthenes(first order reaction)
r(4,1:NC)=k(4)*dpdy(Pa,:)/P;                                                         % P + H2 <=> 2G  [mol/(kg*s)]   Hydrocracking of paraffins (first order reaction)

% Unit change from [kmol/(kg*h)] to [kmol/(kg*s)]
r=r./3600;

return