Hledání optimální teploty pro průběh silně exotermní vratné chemické reakce

%
% OPTIMALIZACE TEPLOTY U SILNE EXOTERMICKE REAKCE
%
% Zadani a predpoklady:
% - vsadkovy reaktor [doba reakce znama]
% - izotermicky prubeh reakce [hledame optimalni teplotu]
% - konstantni objem smesi (reaktanty i produkty rozpustene)
% - vratna reakce aA <---> bB
 
  tenda = 10*60;       % doba reakce A
  tendb = 60*60;       % doba reakce B
 
  cA0 = 1; % mol/dm3
  a   = 1;
  b   = 1;
  R   = 8.314;
 
% kinetika reakce
  A = 2e3;
  Ea  = 40e3;    % aktivacni energie (citlivost na teplote), J/mol
% rovnovaha reakce
  dHr = -70e3;   % reakcni entalpie (exotermicka reakce), J/mol
  K0  = 200;     % rovnovazna konstanta pri T0
  T0  = 298.15;  % K
  Z = exp(log(K0) + dHr/R/T0);
 
% FUNKCE PRO VYPOCET ZAVISLOSTI "K" a "kp" na teplote
  FK  = @(T) Z*exp(-dHr/R/T); % rovnovazna konstanta
  Fkp = @(T) A*exp(-Ea/R/T);  % rychlostni konstanta
 
% FUNKCE PRO VYPOCET ZAVISLOSTI RYCHLOSTI REAKCE na konverzi a teplote
  FR = @(X,T) Fkp(T)*(-cA0^a*(1-X)^a + 1/FK(T)*cA0^b*(b/a*X)^b);
 
% prohledavany rozsah teplot
  T = linspace(T0,450,70);
 
  Xenda = []; % dosazena konverze za dobu A
  Xendb = []; % dosazena konverze za dobu B
  Xeq   = []; % konverze pri rovnovaze
  for Ti=T
    fun = @(t,y) -a*FR(y,Ti)/cA0; % derivace konverze (z molarni bilance)
 
    [tt,yy] = ode15s(fun,[0 tenda], 0);
    Xenda = [Xenda, yy(end)];
    [tt,yy] = ode15s(fun,[0 tendb], 0);
    Xendb = [Xendb, yy(end)];
 
    Xeqnow = fzero(@(X) FR(X,Ti), [0.0000000000001,0.999999999999]);
    Xeq = [Xeq, Xeqnow];
  end
 
  plot(T,Xenda,'o',T,Xendb,'+',T,Xeq,'-')
  xlabel('temperature / K')
  ylabel('conversion / [-]')
  legend('time A','time B','equilibrium')
  title('Batch reactor')