% run_axial_initial2.m
% This program employs the initial value technique
% to solve the reactor problem with axial dispersion.
% The iteration is performed
% using the Newton-Raphson method.
% This goes from the outlet (z = 0) to the inlet (z = 1).
global Pe Da iwhich

% set Pe and Da and iwhich before calling program
tol = 1.e-12;
change = 1.;
al = 0.1;
iter=0;
while change>tol
iter=iter+1;
y0 = [al 0 1 0];
[z c] = ode45('axialrhs2',[0 1],y0);
nn= size(c);
psi = c(nn(1),2)/Pe - 1 + c(nn(1),1);
dpsi = c(nn(1),4)/Pe + c(nn(1),3);
change = abs(psi/dpsi)
al = al - psi/dpsi
end
y0 = [al 0 1 0];
[z c] = ode45('axialrhs2',[0 1],y0);
plot(1-z,c(:,1),'b')
xlabel('x = 1 - z')
ylabel('c')
iter
psi

******************************************************
% axialrhs2.m
% This function evaluates the right-hand side for
% the reactor with axial dispersion.

function ydot=axialrhs(x,y)
global Pe Da iwhich

yyy = y(1);
ans=rate(yyy,iwhich);
rxnrate = ans(1);
drate = ans(2);
ydot(1) = y(2);
ydot(2) = Pe*(-y(2) + Da*rxnrate);
ydot(3) = y(4);
ydot(4) = Pe*(-y(4) + Da*drate*y(3));
ydot = ydot';

******************************************************
% 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