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;
k1 = 1.0e-1; % 1/sec
%  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);
   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);
   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


