function dydt = syspde_para_input(n_eq,i_x_v,y1,t,dxi,mx)
% evaluate first and second spatial derivatives for all variables
for k_eq = 1:1:n_eq
   for i_x = i_x_v(k_eq,1):1:i_x_v(k_eq,2)
%   	dydx_matrix(k_eq, i_x) = 0.5*( y1(k_eq,i_x+1) - y1(k_eq,i_x-1) )*dxi;
%   	dydx_matrix(k_eq, i_x) = ( y1(k_eq,i_x+1) - y1(k_eq,i_x) )*dxi;
   	dydx_matrix(k_eq, i_x) = ( y1(k_eq,i_x) - y1(k_eq,i_x-1) )*dxi;
      d2ydx2_matrix(k_eq,i_x) = ( y1(k_eq,i_x+1) - 2.0*y1(k_eq,i_x) + y1(k_eq,i_x-1) )*dxi^2;
   end
end
%
%  Heterogeneous Catalysis 
%  EXAMPLE 2.  A and B flowing into a pipe packed with adsorbing particles
%              rate constants for adsorption and desorption are known.
%              isothermal, surface adsorption only.
%              Langmuir-Hinshelwood  kinetics 
%              high adsorption capacity relative to bulk density
%              pipe initially filled with solvent
%              diffusion added
%
%
% input parameters
%
%
%  void fraction
eps = 0.50;
%
% pellet diameter
%
Dp = 1.0e-6; % m
%
% site surface density
%
rhoS = 1.0e+18; % sites/m^2
%
%  single spherical particle surface area and volume
%
pi = 2.0*asin(1.0);
Ap = pi*Dp*Dp;
Vp = pi/6.0*Dp*Dp*Dp;
AoVp = Ap/Vp;
%
%  masses
%
MWA = 0.016; % kg/mole
MWB = 0.016; % kg/mole
MWC = 0.032; % kg/mole
MWS = 0.018; % kg/mole
NAV = 6.022e23; % molecules/mole
mA = MWA/NAV; % kg/molecule
mB = MWB/NAV; % kg
mC = MWC/NAV; % kg
mS = MWS/NAV; % kg
%
%  number of sites per particle
%
Ns = Ap*rhoS;
%
%  number of particles per volume of reactor
%
rhoP = (1-eps)/Vp;
%
%  density of sites per volume of reactor
%
rhoR = Ns*rhoP;   
%
% adsorption rate constants
%
kadsA = 1.0e+0;
kadsB = 1.0e-1;
kadsC = 1.0e-2;
kadsS = 1.0e-3;
kdesA = 1.0e-1;
kdesB = 1.0e-2;
kdesC = 1.0e+0;
kdesS = 1.0e+0;
%
% axial velocity
%
vz = 1.0e-2; % m/sec
%
%  bulk diffusivity
%
DA = 1.0e-5; % m^2/sec
DB = DA; % m^2/sec
DC = DA; % m^2/sec
DS = DA; % m^2/sec
%
% reaction stoichiometry
%
nuA = -1.0;
nuB = -1.0;
nuC = 1.0;
nuS = 0.0;
%
%  reaction rate
%
kr = 1.0e-1;
%
%  evaluate temporal derivative for all equations at interior nodes
%
dydt = zeros(n_eq,mx);
%
%  1.  mass balance on total A
%  2.  mass balance on bulk A
%  3.  mass balance on adsorbed A
%  4.  mass balance on total B
%  5.  mass balance on bulk B
%  6.  mass balance on adsorbed B
%  7.  mass balance on total C
%  8.  mass balance on bulk C
%  9.  mass balance on adsorbed C
%  10.  mass balance on total S
%  11.  mass balance on bulk S
%  12.  mass balance on adsorbed S
%
for i_x = i_x_v(k_eq,1):1:i_x_v(k_eq,2)
%
%  convection terms
%
   convA = -vz*dydx_matrix(2, i_x);
   convB = -vz*dydx_matrix(5, i_x);
   convC = -vz*dydx_matrix(8, i_x);
   convS = -vz*dydx_matrix(11, i_x);
%
%  bulk densities
%
   rhoAbulk = y1(2,i_x);
   rhoBbulk = y1(5,i_x);
   rhoCbulk = y1(8,i_x);
   rhoSbulk = y1(11,i_x);
%
%  bulk mass fractions
%
   wAbulk = rhoAbulk/(rhoAbulk+rhoBbulk+rhoCbulk+rhoSbulk);
   wBbulk = rhoBbulk/(rhoAbulk+rhoBbulk+rhoCbulk+rhoSbulk);
   wCbulk = rhoCbulk/(rhoAbulk+rhoBbulk+rhoCbulk+rhoSbulk);
   wSbulk = rhoSbulk/(rhoAbulk+rhoBbulk+rhoCbulk+rhoSbulk);
%
%  bulk diffusivities
%    
%   diffA = DA*d2wAdz2;
%
   denom_back = y1(2,i_x-1) + y1(5,i_x-1) + y1(8,i_x-1) + y1(11,i_x-1);
   wAbulk_back = y1(2,i_x-1)/denom_back;
   wBbulk_back = y1(5,i_x-1)/denom_back;
   wCbulk_back = y1(8,i_x-1)/denom_back;
   wSbulk_back = y1(11,i_x-1)/denom_back;
%
   denom_front = y1(2,i_x+1) + y1(5,i_x+1) + y1(8,i_x+1) + y1(11,i_x+1);
   wAbulk_front = y1(2,i_x+1)/denom_front;
   wBbulk_front = y1(5,i_x+1)/denom_front;
   wCbulk_front = y1(8,i_x+1)/denom_front;
   wSbulk_front = y1(11,i_x+1)/denom_front;
%
   d2wAdz2 = ( wAbulk_front - 2.0*wAbulk + wAbulk_back )*dxi^2;
   d2wBdz2 = ( wBbulk_front - 2.0*wBbulk + wBbulk_back )*dxi^2;
   d2wCdz2 = ( wCbulk_front - 2.0*wCbulk + wCbulk_back )*dxi^2;
   d2wSdz2 = ( wSbulk_front - 2.0*wSbulk + wSbulk_back )*dxi^2;
%
   diffA = DA*d2wAdz2;
   diffB = DB*d2wBdz2;
   diffC = DC*d2wCdz2;
   diffS = DS*d2wSdz2;
%
%  fractional occupancies
%
   thetaA = y1(3,i_x);
   thetaB = y1(6,i_x);
   thetaC = y1(9,i_x);
   thetaS = y1(12,i_x);
   thetaEmpty = 1.0 - thetaA - thetaB - thetaC - thetaS;
%
%  rates of adsorption and desorption
%
   rateadsA = kadsA*rhoAbulk*thetaEmpty;
   rateadsB = kadsB*rhoBbulk*thetaEmpty;
   rateadsC = kadsC*rhoCbulk*thetaEmpty;
   rateadsS = kadsS*rhoSbulk*thetaEmpty;
   ratedesA = kdesA*thetaA;
   ratedesB = kdesB*thetaB;
   ratedesC = kdesC*thetaC;
   ratedesS = kdesS*thetaS;
%
%  reaction rate
%
   raterxn = kr*thetaA*thetaB;
%
% total balances on A and B
%
   dydt(1,i_x)  = eps*convA + eps*diffA + nuA*raterxn;
   dydt(4,i_x)  = eps*convB + eps*diffB + nuB*raterxn;
   dydt(7,i_x)  = eps*convC + eps*diffC + nuC*raterxn;
   dydt(10,i_x) = eps*convS + eps*diffS + nuS*raterxn;
%
%  bulk balances
%
   dydt(2,i_x) = convA + diffA + 1.0/eps*(-rateadsA + ratedesA);
   dydt(5,i_x) = convB + diffB + 1.0/eps*(-rateadsB + ratedesB);
   dydt(8,i_x) = convC + diffC + 1.0/eps*(-rateadsC + ratedesC);
   dydt(11,i_x) = convS + diffS + 1.0/eps*(-rateadsS + ratedesS);
%
%  adsorbed balances
%
   dydt(3,i_x) = 1.0/(rhoR*mA)*(rateadsA - ratedesA + nuA*raterxn) ;
   dydt(6,i_x) = 1.0/(rhoR*mB)*(rateadsB - ratedesB + nuB*raterxn);
   dydt(9,i_x) = 1.0/(rhoR*mC)*(rateadsC - ratedesC + nuC*raterxn);
   dydt(12,i_x) = 1.0/(rhoR*mS)*(rateadsS - ratedesS + nuS*raterxn);
   if (i_x == 3)
%   fprintf(' theta A B C S %e %e %e %e\n', thetaA, thetaB, thetaC, thetaS); 
%   fprintf(' rhobulk A B C S %e %e %e %e\n', rhoAbulk, rhoBbulk, rhoCbulk, rhoSbulk); 
%   fprintf(' rhotot A B C S %e %e %e %e\n', y1(1,i_x), y1(4,i_x), y1(7,i_x), y1(10,i_x)); 
%   fprintf(' convC diffC rxnC %e %e %e \n', eps*convC, eps*diffC, rhoR*mC*nuC*raterxn ); 
 
  end
end


