clc; clear;
y0(1) = 0.06; % SO2
y0(2) = 0.15; % O2
y0(3) = 0.00; % SO3
y0(4) = 0.79; % N2
F0    = 1.5e-6 ; % Feed rate [mol/s]
key = 1;         % index of the key compound
kr = 1.05e7;  % Preexponential factor [mol3/m3/s]
EA = 60000;   % Activation energy
A1 = -11.12;  % Constant
A2 = 11750;   % Constant
ni = [ -1 -0.5 1 0]; % Stoichiometric matrix
y = @(X)((y0 - ni./ni(key).*y0(key).*X)./(1-1./2.*y0(key).*X));
Tr= [700:10:1050];
Xr = [0 : 0.1 : 1];
[~,nt] = size(Tr);
[~,nr] = size(Xr);
for i = 1 : nr
    yy = y(Xr(i));
    rv_old = 0;
    label{i} = num2str(Xr(i),'X=%3.1f\n');
    for  j = 1:nt
        T = Tr(j);
        k = kr .* exp (-EA./(8.314*T));
        K = exp(A1 + A2./T);
        rv = k .*(yy(1).*sqrt(yy(2)) - yy(3) ./ K);
        result(j,i) = rv;
        if(rv > rv_old)
            maxR(i,1) = T;
            maxR(i,2) = rv;
        end
        rv_old = rv;
    end
end
label{nr+1} = 'maxR';
plot(Tr,result,maxR(:,1),maxR(:,2),'o')
title('Maximum reaction rate in SO_2 oxidation (optimal temp. profile)');
legend(label,'location','NorthEastOutside')
yrange=[-100 300];
xrange=[700 1050];
axis([xrange yrange])
xlabel('T / K')
ylabel('r_V / (mol s^{-1} m^{-3})')
grid on
print -dpng ./maxrate_SO2.png