function IR5cIsak
[x,y]=ode23s(@batch,[0:1:140],[3.996 0 0])
plot(x,y)
eta=(3.996-y(:,1))/3.996;
figure(1);plot(x,y);title('C')
legend('Cxylos','CH2','Cxylitol')
figure(2);plot(x,eta);grid on;title('omsättningsgrad')
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=57.5*1.5;
pB=75000/5200;
ctot=3.997;
p=40;
mcat=75000;  % kg katalysator
T=383.15;
R=0.0831451;
V=5200;

% 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;
vatgasflode=(cG-Keq*cB)/(Keq/kLa)*V

return