function [f,F,diff,x_var,p_var,d_var,L] = CSTR_model_ode(par)
% Written by Krishnamoorthy, Oct 2018.
% email - dineshk@ntnu.no
%
% Original model adapted from Srinivasan and Bonvin, 2007. 
%
import casadi.*

CA = MX.sym('CA',1);
CB = MX.sym('CB',1);
CC = MX.sym('CC',1);
qr = MX.sym('qr',1);

nA = MX.sym('nA',1);
nB = MX.sym('nB',1);
nC = MX.sym('nC',1);

FA = MX.sym('FA',1);
FB = MX.sym('FB',1);

k1 = MX.sym('k1',1);

CA = nA/(par.V);
CB = nB/(par.V);
CC = nC/(par.V);

dx1 = (FA*par.CAin - (FA+FB)*CA - k1*CA*CB*par.V);
dx2 = (FB*par.CBin - (FA+FB)*CB - k1*CA*CB*par.V - 2*par.k2*CB^2*par.V);
dx3 = (-(FA+FB)*CC + k1*CA*CB*par.V);

diff = vertcat(dx1,dx2,dx3);

x_var = vertcat(nA,nB,nC);
p_var = vertcat(FA,FB);
d_var = k1;

L = -((FA+FB)^2*CC^2)/(FA*par.CAin);

ode = struct('x',x_var,'p',vertcat(p_var,k1),'ode',diff,'quad',L);  
opts = struct('tf',par.Ts/par.t);

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

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