function syspde_para
clear all;
close all;
% n_eq = number of equations
n_eq = 12; %INPUT

% m_int = number of spatial intervals
m_int = 20; %INPUT
xo = 0.0;
xf = 1.0;
n_int = 2000;
to = 0.0;
tf = 200.0;
%  BC is a vector of boundary condition type definitions
%  0 = Dirichlet
%  1 = Neumann
BC = [0 1; 0 1; 0 1; 0 1; 0 1; 0 1; 0 1; 0 1; 0 1; 0 1; 0 1; 0 1]; %INPUT

%  code must create files
% function y_out = syspde_boundary_conditions(k_eq,k_bc,x,t,y1,i_x, dx);
% function y_out = syspde_initial_conditions(k_eq,x);





%
%  create spatial grid
%
%  mx = number of spatial gridpoints
%  (this grid has imaginary points for Neumann BCs)
%  (it still works if the problem has only Dirichlet BCs)
mx = m_int + 3;
dx = (xf-xo)/m_int;
dxi = 1/dx;
xgrid = [(xo-dx):dx:(xf+dx)]';
%  create temporal grid
%  nt = number of temporal gridpoints
nt = n_int + 1;
dt = (tf-to)/n_int;
tgrid = [(to):dt:(tf)]';
%
%  dimension solution matrix
%
y = zeros(n_eq,nt,mx);
%
%  define some limitations on indices for where solution is not given by ICs and
%  where the derivatives need to be calculated using finite difference techniques
%
for k_eq = 1:1:n_eq
   for i_BC = 1:1:2
      if (i_BC == 1)
   		if (BC(k_eq, i_BC) == 0) 
      		i_x_v(k_eq,i_BC) = 3;
   		else
      		i_x_v(k_eq,i_BC) = 2;
         end
      else
   		if (BC(k_eq, i_BC) == 0) 
      		i_x_v(k_eq,i_BC) = mx-2;
   		else
      		i_x_v(k_eq,i_BC) = mx-1;
         end
      end
   end
end
%i_x_v
%i_x_d
%
%  input initial conditions into solution matrix
%
j_t = 1;
for k_eq = 1:1:n_eq
   for i_x = 2:1:mx-1
      x = xgrid(i_x);
      y(k_eq,j_t,i_x) = syspde_initial_conditions(k_eq,x);
   end
end
%
%  input Dirichlet Boundary conditions into solution matrix for all times
%
y1 = zeros(n_eq,mx);
for k_eq = 1:1:n_eq
   for k_bc = 1:1:2
   	if (BC(k_eq,k_bc) == 0) 
         if (k_bc == 1) 
            i_x = 2;
         else
            i_x = mx - 1;
         end
			x = xgrid(i_x);
         for j_t = 1:1:nt
            t = tgrid(j_t);
            for m_eq = 1:1:n_eq
               y1(m_eq,:) = y(m_eq,j_t,:);
            end
%            fprintf(1,'D: k_eq= %i kbc= %i i_x = %i x = %f t = %f  \n',k_eq, k_bc,i_x,x,t)
            y(k_eq,j_t,i_x) = syspde_boundary_conditions(k_eq,k_bc,x,t,y1,i_x,dx);
         end
      end
   end
end
%
%  input initial Neumann Boundary conditions into solution matrix
%
j_t = 1;
t = tgrid(j_t);
for k_eq = 1:1:n_eq
	y1(k_eq,:) = y(k_eq,j_t,:);
end
for k_eq = 1:1:n_eq
   for k_bc = 1:1:2
   	if (BC(k_eq,k_bc) == 1) 
         if (k_bc == 1) 
            i_x = 1;
            i_x2 = 2;
            i_x3 = 3;
         else
            i_x = mx;
            i_x2 = mx-1;
            i_x3 = mx-2;
         end
         x = xgrid(i_x);
%         fprintf(1,'N: k_eq= %i kbc= %i i_x = %i x = %f t = %f  \n',k_eq, k_bc,i_x,x,t)
			y(k_eq,j_t,i_x) = syspde_boundary_conditions(k_eq,k_bc,x,t,y1,i_x,dx);
      end
   end
end
%
%  call Classical second order Runge-Kutta method 
%
dydt_1 = zeros(n_eq,mx);
dydt_2 = zeros(n_eq,mx);
for j_t = 2:1:nt
   %  obtain first half of approximation to derivative 
   t = tgrid(j_t-1);
   for k_eq = 1:1:n_eq
		y1(k_eq,:) = y(k_eq,j_t-1,:);
	end
   dydt_1 = syspde_para_input(n_eq,i_x_v,y1,t,dxi,mx);
   %  obtain second half of approximation to derivative 
   t = t + dt;
      %  first obtain intermediate values of y in interior nodes
   for k_eq = 1:1:n_eq
      for i_x = i_x_v(k_eq,1):1:i_x_v(k_eq,2)
         y1(k_eq,i_x) = y1(k_eq,i_x) + dt*dydt_1(k_eq,i_x);
      end
   end
   %  second obtain intermediate values of y from Dirichlet BCs
   for k_eq = 1:1:n_eq
   	for k_bc = 1:1:2
   		if (BC(k_eq,k_bc) == 0) 
          	if (k_bc == 1) 
            	i_x = 2;
         	else
            	i_x = mx-1;
	         end
         	x = xgrid(i_x);
            y1(k_eq,i_x) = y(k_eq,j_t,i_x);
         end
      end
   end
   % third obtain intermediate values of y from Neumann BC's
	t = tgrid(j_t);
	for k_eq = 1:1:n_eq
   	for k_bc = 1:1:2
   		if (BC(k_eq,k_bc) == 1) 
         	if (k_bc == 1) 
            	i_x = 1;
         	else
            	i_x = mx;
	         end
         	x = xgrid(i_x);
				y1(k_eq,i_x) = syspde_boundary_conditions(k_eq,k_bc,x,t,y1,i_x,dx);
      	end
   	end
	end
   dydt_2 = syspde_para_input(n_eq,i_x_v,y1,t,dxi,mx);
   for k_eq = 1:1:n_eq
      for i_x = 1:1:mx
         y(k_eq,j_t,i_x) = y(k_eq,j_t-1,i_x) + 0.5*dt*(dydt_1(k_eq,i_x) + dydt_2(k_eq,i_x) );
      end
   end
   % compute new values of unknown from Neumann BC's
	t = tgrid(j_t);
	for k_eq = 1:1:n_eq
   	for k_bc = 1:1:2
   		if (BC(k_eq,k_bc) == 1) 
         	if (k_bc == 1) 
            	i_x = 1;
         	else
            	i_x = mx;
	         end
         	x = xgrid(i_x);
				y(k_eq,j_t,i_x) = syspde_boundary_conditions(k_eq,k_bc,x,t,y1,i_x,dx);
      	end
   	end
	end
end
%
%  write solution
%
%y
%
%  plot solution
%
% define plot/movie parameters
T=y;
nframe = 1000;
if (n_int > nframe)
	nskip = round(n_int/nframe);
else
   nskip = 1;
end
lplot = 1;
fps = 3;
%  define plot line colors (for up to 10 eqns)
plc{1} = 'k-';
plc{2} = 'k:';
plc{3} = 'k--';
plc{4} = 'r-';
plc{5} = 'r:';
plc{6} = 'r--';
plc{7} = 'b-';
plc{8} = 'b:';
plc{9} = 'b--';
plc{10} = 'g-';
plc{11} = 'g:';
plc{12} = 'g--';
plc{13} = 'm-';
plc{14} = 'm:';
plc{15} = 'm--';
if (lplot == 1)
   fpsinv = 1.0/fps;
	newplot;
	Tmax=max(max(max(T)));
   Tmin=min(min(min(T)));
   xtext = xo + 0.1*(xf-xo);
   ytext = Tmin +0.9*(Tmax-Tmin);
   for j_t=1:nskip:nt
      for k_eq = 1:1:n_eq
         for i_x = 2:1:mx-1
            vector_plot(i_x) = T(k_eq,j_t,i_x);
         end
         plot(xgrid(2:mx-1) , vector_plot(2:mx-1),plc{k_eq} );
         hold on
      end
      axis([xo xf Tmin Tmax]);
      XLABEL('position along rod (cm)') 
      YLABEL('Temperature (K)') 
      temps = char(strcat('time = ', num2str(tgrid(j_t)), ' sec '));
      text(xtext, ytext,temps)
      hold off
      pause(fpsinv);
	end
end

function y_out = syspde_boundary_conditions(k_eq,k_bc,x,t,y1,i_x, dx);
%
%  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;   
Ctot = 40.6; % moles/m^3
xA = 0.95;
xB = 0.05;
xC = 0.0;
xS = 1.0- xA - xB - xC;
%
%  bulk densities
%
rhobulkA = Ctot*xA*MWA; % kg/m^3
rhobulkB = Ctot*xB*MWB; % kg/m^3
rhobulkC = Ctot*xC*MWC; % kg/m^3
rhobulkS = Ctot*xS*MWS; % kg/m^3
%
% adsorption rate constants
%
kadsA = 1.0e-2;
kadsB = 1.0e+0;
kadsC = 1.0e-2;
kadsS = 1.0e-3;
kdesA = 1.0e+0;
kdesB = 1.0e-2;
kdesC = 1.0e+0;
kdesS = 1.0e+0;
%
%  adsorption equilibrium coefficients
%
KA = kadsA/kdesA;
KB = kadsB/kdesB;
KC = kadsC/kdesC;
KS = kadsS/kdesS;
%
%  fractional occupancies
%
thetaA = KA*rhobulkA/(1.0 + KA*rhobulkA + KB*rhobulkB + KC*rhobulkC + KS*rhobulkS);
thetaB = KB*rhobulkB/(1.0 + KA*rhobulkA + KB*rhobulkB + KC*rhobulkC + KS*rhobulkS);
thetaC = KC*rhobulkC/(1.0 + KA*rhobulkA + KB*rhobulkB + KC*rhobulkC + KS*rhobulkS);
thetaS = KS*rhobulkS/(1.0 + KA*rhobulkA + KB*rhobulkB + KC*rhobulkC + KS*rhobulkS);
%
%  total density
%
rhototA = eps*rhobulkA + mA*rhoR*thetaA;
rhototB = eps*rhobulkB + mB*rhoR*thetaB;
rhototC = eps*rhobulkC + mC*rhoR*thetaC;
rhototS = eps*rhobulkS + mS*rhoR*thetaS;
%
%  specify boundary conditions
%
if (k_bc == 1)
	if (k_eq == 1) 
   	      y_out = rhototA;
	elseif (k_eq == 2)
		y_out = rhobulkA;
	elseif (k_eq == 3)
		y_out = thetaA;
	elseif (k_eq == 4)
   	      y_out = rhototB;
	elseif (k_eq == 5)
		y_out = rhobulkB;
	elseif (k_eq == 6)
		y_out = thetaB;
	elseif (k_eq == 7)
   	      y_out = rhototC;
	elseif (k_eq == 8)
		y_out = rhobulkC;
	elseif (k_eq == 9)
		y_out = thetaC;
	elseif (k_eq == 10)
   	      y_out = rhototS;
	elseif (k_eq == 11)
		y_out = rhobulkS;
	elseif (k_eq == 12)
		y_out = thetaS;
   end
else
	if (k_eq == 1) 
		y_out = y1(k_eq, i_x-2)- 2*dx*(0.0);
	elseif (k_eq == 2)
		y_out = y1(k_eq, i_x-2)- 2*dx*(0.0);
	elseif (k_eq == 3)
		y_out = y1(k_eq, i_x-2)- 2*dx*(0.0);
	else
		y_out = y1(k_eq, i_x-2)- 2*dx*(0.0);
   end
end
%fprintf(1,'BC: k_eq= %i kbc= %i i_x = %i x = %f t = %f yout = %f \n',k_eq, k_bc,i_x,x,t,y_out)

function y_out = syspde_initial_conditions(k_eq,x);
%
%  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;   
Ctot = 40.6; % moles/m^3
xA = 0.0;
xB = 0.0;
xC = 0.0;
xS = 1.0- xA - xB - xC;
%
%  bulk densities
%
rhobulkA = Ctot*xA*MWA; % kg/m^3
rhobulkB = Ctot*xB*MWB; % kg/m^3
rhobulkC = Ctot*xC*MWC; % kg/m^3
rhobulkS = Ctot*xS*MWS; % kg/m^3
%
% adsorption rate constants
%
kadsA = 1.0e-2;
kadsB = 1.0e+0;
kadsC = 1.0e-2;
kadsS = 1.0e-3;
kdesA = 1.0e+0;
kdesB = 1.0e-2;
kdesC = 1.0e+0;
kdesS = 1.0e+0;
%
%  adsorption equilibrium coefficients
%
KA = kadsA/kdesA;
KB = kadsB/kdesB;
KC = kadsC/kdesC;
KS = kadsS/kdesS;
%
%  fractional occupancies
%
thetaA = KA*rhobulkA/(1.0 + KA*rhobulkA + KB*rhobulkB + KC*rhobulkC + KS*rhobulkS);
thetaB = KB*rhobulkB/(1.0 + KA*rhobulkA + KB*rhobulkB + KC*rhobulkC + KS*rhobulkS);
thetaC = KC*rhobulkC/(1.0 + KA*rhobulkA + KB*rhobulkB + KC*rhobulkC + KS*rhobulkS);
thetaS = KS*rhobulkS/(1.0 + KA*rhobulkA + KB*rhobulkB + KC*rhobulkC + KS*rhobulkS);
%
%  total density
%
rhototA = eps*rhobulkA + mA*rhoR*thetaA;
rhototB = eps*rhobulkB + mB*rhoR*thetaB;
rhototC = eps*rhobulkC + mC*rhoR*thetaC;
rhototS = eps*rhobulkS + mS*rhoR*thetaS;
%
%  specify initial conditions
%
	if (k_eq == 1) 
   	      y_out = rhototA;
	elseif (k_eq == 2)
		y_out = rhobulkA;
	elseif (k_eq == 3)
		y_out = thetaA;
	elseif (k_eq == 4)
   	      y_out = rhototB;
	elseif (k_eq == 5)
		y_out = rhobulkB;
	elseif (k_eq == 6)
		y_out = thetaB;
	elseif (k_eq == 7)
   	      y_out = rhototC;
	elseif (k_eq == 8)
		y_out = rhobulkC;
	elseif (k_eq == 9)
		y_out = thetaC;
	elseif (k_eq == 10)
   	      y_out = rhototS;
	elseif (k_eq == 11)
		y_out = rhobulkS;
	elseif (k_eq == 12)
		y_out = thetaS;
   end