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
%
%  EXAMPLE 2.  A  flowing into a pipe initially filled with Solvent, S.
%  elementary reaction A --> B
%  Isothermal.
%
%
% input parameters
%
% axial velocity
%
vz = 1.0e-2; % m/sec
%
%  diffusivity
%
DA = 1.0e-9; % m^2/sec
DS = DA; % m^2/sec
DB = DA; % m^2/sec
%
% stoichiometry
%
nuA = -1;
nuB = 1;
ko1 = 1.0e-1; % 1/sec
R = 8.314; % J/mol/K
Ea1 = 2000; % J/mol
%
%  heat capacities
%
Tref = 298.15; % K
CpA = 80.0; % J/mol/K
CpB = 100.0; % J/mol/K
CpS = 50.0; % J/mol/K
%
%  thermal conductivity
%
kc = 6.0e+0; % J/s/m/K
%
%  heat of reaction
%
dHR = -50000.0; % J/mol
%
%  heat loss
%
% diameter
Dp = 0.1; % m
% surface area per length
Asurf = pi*Dp;
% reactor volume per length
Vr = 0.25*pi*Dp*Dp;
% heat transfer coefficient
h = 0.0e0; % J/s/m^2/K
%  surrounding temperature
Tsurr = 300.0; % K

%  evaluate temporal derivative for all equations at interior nodes
dydt = zeros(n_eq,mx);
%
%  mole balance on A
%
k_eq = 1;
for i_x = i_x_v(k_eq,1):1:i_x_v(k_eq,2)
   conv = -vz*dydx_matrix(k_eq, i_x);
   diff = + DA*d2ydx2_matrix(k_eq,i_x);
   CA = y1(1,i_x);
   T  = y1(4,i_x);
   k1 = ko1*exp(-Ea1/(R*T));
   rate = k1*CA^abs(nuA);
   dydt(k_eq,i_x) = conv + diff + nuA*rate;
end
%
%  mole balance on B
%
k_eq = 2;
for i_x = i_x_v(k_eq,1):1:i_x_v(k_eq,2)
   conv = -vz*dydx_matrix(k_eq, i_x);
   diff = + DA*d2ydx2_matrix(k_eq,i_x);
   CA = y1(1,i_x);
   T  = y1(4,i_x);
   k1 = ko1*exp(-Ea1/(R*T));
   rate = k1*CA^abs(nuA);
   dydt(k_eq,i_x) = conv + diff + nuB*rate;
end

%
%  mole balance on S
%
k_eq = 3;
for i_x = i_x_v(k_eq,1):1:i_x_v(k_eq,2)
   conv = -vz*dydx_matrix(k_eq, i_x);
   diff = + DS*d2ydx2_matrix(k_eq,i_x);
   dydt(k_eq,i_x) = conv + diff;
end

%
%  energy balance 
%
k_eq = 4;
for i_x = i_x_v(k_eq,1):1:i_x_v(k_eq,2)
   CA = y1(1,i_x);
   CB = y1(2,i_x);
   CS = y1(3,i_x);
   T  = y1(4,i_x);
   dCAdz = dydx_matrix(1, i_x);
   dCBdz = dydx_matrix(2, i_x);
   dCSdz = dydx_matrix(3, i_x);
   dTdz  = dydx_matrix(4, i_x);
   dCAdt = dydt(1, i_x);
   dCBdt = dydt(2, i_x);
   dCSdt = dydt(3, i_x);
   k1 = ko1*exp(-Ea1/(R*T));
   rate = k1*CA^abs(nuA);
   sumCpdCdz = CpA*dCAdz + CpB*dCBdz + CpS*dCSdz;
   sumCpC = CpA*CA + CpB*CB + CpS*CS;
   sumCpdCdt = CpA*dCAdt + CpB*dCBdt + CpS*dCSdt;
   term1 = -vz*(T-Tref)*sumCpdCdz;
   term2 = -vz*dTdz*sumCpC;
   term3 = kc*dydx_matrix(k_eq, i_x);
   term4 = -dHR*rate;
   term5 = Asurf*h*(Tsurr-T)/Vr;
   term6 = -(T-Tref)*sumCpdCdt;
   EB_numerator = term1 + term2 + term3 + term4 + term5 + term6;
   EB_denominator = sumCpC;
   dydt(k_eq,i_x) = EB_numerator/EB_denominator;
end

