% filename rxnFD.m
% This code solves the equations for reaction-diffusion
% in a porous catalyst using the finite difference method.
% a = 1, 2, 3 for planar, cylindrical, and spherical geometry
% If the problem is linear, set linear = 1 to avoid a second iteration.

geometry = 3
linear = 0
%set ngrid before calling
% ngrid must be a number*2+1 (for Simpson's rule)
% set phi before calling %phi = 10
phi2 = phi*phi;
dr = 1/(ngrid-1);
dr2=dr*dr;

% iterate
tol = 1.e-12;
change = 1.;
iter=0;
% initial guess
for i=1:ngrid
c(i) = 0.0;
rad(i) = (i-1)/(ngrid-1);
end

while change>tol
iter=iter+1;

% Set up the matrices
% set the matrices to zero
for i=1:ngrid
f(i) = 0.;
for j=1:ngrid
aa(i,j) = 0.;
end
end
% set the matrices for the differential equation
for j=2:ngrid-1
ans = rate(c(j));
raterxn = ans(1);
drate = ans(2);
aa(j,j-1) = 1 - 2.*dr/rad(j);
aa(j,j) = -2 - phi2*dr2*drate;
aa(j,j+1) = 1 + 2.*dr/rad(j);
f(j) = phi2*dr2*(raterxn - drate*c(j));
end
f(ncol+1) = 1;
ans = rate(c(1));
aa(1,1) = -1 - dr2*phi2*ans(2)/6.;
aa(1,2) = 1;
f(1) = phi2*dr2*(ans(1) - ans(2)*c(1))/6.;
aa(ngrid,ngrid) = 1.;
f(ngrid)=1.;

% Solve one iteration.
c11 = aa\f';

% Calculate the criterion to stop the iterations.
sum=0.;
for i=1:ngrid
sum = sum + abs(c11(i) - c(i));
c(i) = c11(i);
end
if linear==1,break,end
change = sum
end

% Plot the solution.
plot(rad,c,'*-b')
xlabel('r')
ylabel('c')
c

% calculate the effectiveness factor
for i=2:2:ngrid-1
w(i) = 4.;
w(i+1) = 2.;
end
w(1) = 1;
w(ngrid) = 1;
sum = 0.;
for i=1:ngrid
ans = rate(c(i));
raterxn = ans(1);
num1 = rad(i)^(geometry-1);
sum = sum + w(i)*raterxn*num1;
end
eta = geometry*(dr/3)*sum



************************************************
% rate.m
% calculate the rate and derivative for a given concentration

function ans = rate(c)
% use for linear reaction
%ans(1) = c;
%ans(2) = 1.;

% use for second order reaction
%ans(1) = c*c;
%ans(2) = 2.*c;

% use for non-isothermal reaction
%ggg = 18.; bbb = 0.3;
ggg = 30; bbb = 0.4;
t = 1 + bbb - bbb*c;
eee = exp(ggg*(1 - 1/t));
ans(1) = c*eee;
ans(2) = eee*(1 - bbb*ggg*c/(t*t));