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

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

par.Ts = 10;

dist = [1;0.0];

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

par.lbu = [0;0;0];
par.ubu = [1000;2;2];

par.dx0 = [0.4978;0.5022;424];
par.u0 = [420.2;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);

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

E1 = MX.sym('E1',1);
E2 = MX.sym('E2',1);

du = 1;
% ODE system
dx1 = (CAin - CA)/tau - (A1*exp(-E1/(R*T))*CA - A2*exp(-E2/(R*T))*CB);
dx2 = (CBin - CB)/tau + (A1*exp(-E1/(R*T))*CA - A2*exp(-E2/(R*T))*CB);
dx3 = (Tin - T)/tau + 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+var*par.E2);

x_var = vertcat(CA,CB,T);
p_var = vertcat(Tin,CAin,CBin);
d_var = vertcat(CAin,CBin);
L = -(pCB*CB - (pTi*Tin).^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 = 1;


%% Steady state optimization
dist = [1;0];
w = {};
w0 = [];
lbw = [];
ubw = [];

g = {};
lbg = [];
ubg = [];

w = {w{:},x_var,p_var};
lbw = [lbw;par.lbx;par.lbu(1);dist];
ubw = [ubw;par.ubx;par.ubu(1);dist];
w0 = [w0;par.dx0;par.u0(1);dist];

J = L; % Economic objective

g = {g{:},vertcat(diff)};
lbg = [0.*ones(3,1)];
ubg = [0.*ones(3,1)];

nlp = struct('x',vertcat(w{:}),'f',J,'g',vertcat(g{:}));
solver = nlpsol('solver','ipopt',nlp);
sol = solver('x0',w0,'lbx',lbw,'ubx',ubw,'lbg',lbg,'ubg',ubg);
w_opt_SS = full(sol.x);
disp('-------Steady state Optimum-------')
disp(' ')
display([w_opt_SS(1:4);-full(sol.f)]);
disp('----------------------------------')

%%
diffT = (Hr*((A1*exp(-E1/(R*T)))*CA - (A2*exp(-E2/(R*T)))*CB)/(rho*cp));
diffT = substitute(diffT,A1,par.A1);
diffT = substitute(diffT,A2,par.A2);
diffT = substitute(diffT,E1,par.E1);
diffT = substitute(diffT,E2,par.E2);
df_dT= Function('dJ_du',{x_var,p_var},{jacobian(diffT,x_var(3))});