% filename axialFD.m
% This code solves the equations for a reactor with
% axial dispersion 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.

% set linear to 1 if problem is linear - this avoids one unnecessary iteration
linear = 0
% set the type of reaction rate
iwhich = 2
%set ngrid, Peclet (Pe), and Damkohler (Da) before calling
peinv = 1/Pe;
dx = 1/(ngrid-1);
dx2=dx*dx;

% iterate
tol = 1.e-12;
change = 1.;
iter=0;
% initial guess
for i=1:ngrid
c(i) = 0.0;
x(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),iwhich);
raterxn = ans(1);
drate = ans(2);
aa(j,j-1) = peinv + dx/2;
aa(j,j) = -2*peinv - Da*dx2*drate;
aa(j,j+1) = peinv - dx/2;
f(j) = Da*dx2*(raterxn - drate*c(j));
end
ans = rate(c(1),iwhich);
raterxn = ans(1);
drate = ans(2);
aa(1,1) = -2*peinv - 2*dx - Pe*dx2 - Da*dx2*drate;
aa(1,2) = 2*peinv;
f(1) = Da*dx2*(raterxn - drate*c(1)) - 2*dx - dx2*Pe;
ans = rate(c(ngrid),iwhich);
raterxn = ans(1);
drate = ans(2);
aa(ngrid,ngrid-1) = 2*peinv;
aa(ngrid,ngrid) = -2*peinv -Da*dx2*drate;
f(ngrid)=Da*dx2*(raterxn - drate*c(ngrid));

% 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(x,c,'*-b')
axis([0 1 0 1])
xlabel('r')
ylabel('c')
c

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

function ans = rate(c,iwhich)

if iwhich==1
% use for linear reaction
ans(1) = c;
ans(2) = 1.;

elseif iwhich == 2
% use for second order reaction
ans(1) = c*c;
ans(2) = 2.*c;
else
% 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));
end