function [par,diff,x_var,u_var,p_var,d_var,L,F_integrator,f] = CSTR_model(par)
import casadi.*

par.A1  = 5000;
par.A2  = 1e6;
cp      = 1000;
par.E1  = 10000;
par.E2  = 15000;
R       = 1.987;
Hr      = 5000;
rho     = 1;
V       = 1*60;
pCB     = 2.009;
pTi     = 1.657e-3;
F_sp    = 1;

par.Ts  = 1;

par.lbx     = [0;0;0];
par.ubx     = [100;100;1000];

par.lbu     = [0;1;par.dist];
par.ubu     = [1000;1;par.dist];

par.dx0     = [0.4978;0.5022;426.71];
par.u0      = [424.2;1;par.dist];

% x
CA      = MX.sym('CA',1);
CB      = MX.sym('CB',1);
T       = MX.sym('T',1);

% dist
CAin    = MX.sym('CAin',1);
CBin    = MX.sym('CBin',1);

% u
Tin     = MX.sym('Tin',1);
F       = MX.sym('F',1);

% p
A1      = MX.sym('A1',1);
A2      = MX.sym('A2',1);
E1      = MX.sym('E1',1);
E2      = MX.sym('E2',1);

% ODE system
dx1     = F*(CAin - CA)/V - (A1*exp(-E1/(R*T))*CA - A2*exp(-E2/(R*T))*CB);
dx2     = F*(CBin - CB)/V + (A1*exp(-E1/(R*T))*CA - A2*exp(-E2/(R*T))*CB);
dx3     = F*(Tin - T)/V + Hr*((A1*exp(-E1/(R*T)))*CA - (A2*exp(-E2/(R*T)))*CB)/(rho*cp);

diff    = vertcat(dx1,dx2,dx3);

diff    = substitute(diff,A1,par.A1);
diff    = substitute(diff,A2,par.A2);
diff    = substitute(diff,E1,par.E1);
diff    = substitute(diff,E2,par.E2);

x_var   = vertcat(CA,CB,T);
u_var   = vertcat(Tin,F);
p_var   = vertcat(u_var,CAin,CBin);
d_var   = vertcat(CAin,CBin);
L       = -(pCB*CB - (pTi*Tin).^2);% + 10000*(F - F_sp)^2;

ode     = struct('x',x_var,'p',p_var,'ode',diff,'quad',L);  % der(m_tot) = w_in - w_out;
opts    = struct('tf',par.Ts);

f       = Function('f',{x_var,p_var},{diff,L},{'x','p'},{'xdot','qj'});

% create IDAS integrator
F_integrator = integrator('F','cvodes',ode,opts);

par.nx  = numel(x_var);
par.nd  = numel(d_var);
par.nu  = numel(u_var);


