function dx = model(t,x)

%Parameter
K = 1;
kAB0 = 0.5;
kAB = kAB0 * (1-exp(-K*x(1)));% + ;
kBC=3;
kCB=2;
Vdot=1;
VR=10;
ca0 = 1;

dx = zeros(3,1);

dx(1) = Vdot/VR*(ca0 - x(1)) - kAB*x(1);
dx(2) = -Vdot/VR * x(2) + kAB*x(1) + kCB*x(3) - kBC*x(2)^2;
dx(3) = -Vdot/VR * x(3) - kCB*x(3) + kBC*x(2)^2;



