function IR5Isak
[x,y]=ode23s(@batch,[0 120],[4.307 0 0])
plot(x,y)
eta=(4.307-y(:,1))/4.307
figure(1);plot(x,y)
figure(2);plot(x,eta);grid on
figure(3);plot(x,y(:,2));title('vatgas koncentration')
return

function dcdt = batch(t,c)

% rate equations
k  = 2.8;
Kxyl  = 5.26e-3;
KH2  = 0;


cA = c(1);
cB = c(2);
cC = c(3);


%kLa=0.1e10; % intensiv omr?rning
kLa=1000000;
pB=75000/5200;
ctot=4.3089;
p=40;
mcat=75000;  % kg katalysator
T=383.15;
R=0.0831451;

% gas-liquid equilibrium constants
cG=p/(R*T);

xH2=(cA/(cA+cC))*(log(0.9991)-(0.1144/T)+0.0004228*log(p))+(cC/(cA+cC))*(log(0.9993)-(0.1603/T)+0.00041126*log(p));
Keq = cG/(xH2*ctot);
Rate=(k*cA*cB)/((1+Kxyl*cA+sqrt(KH2*cB))^3);

% mass balances
dcdt(1,1) = -pB*Rate;
dcdt(2,1) = (cG-Keq*cB)/(Keq/kLa)-0.8*pB*Rate;
dcdt(3,1) = pB*Rate;

return