function [timeElapsed, uOpt, obj, P, M, iterations, success, exitFlag, x1Horizon, firstx2, lastx2] = MultipleShootingLabview(x0, u0, Qin, Betain, h, Q, R, x2SP, Qu, u_min, u_max, du_max, mh, ph, deltaX, P2, valveConstant)

tic % Starting timer
import casadi.*

%% Model variables

x1 = SX.sym('x1');   % Overflow oil fraction, BetaO
x2 = SX.sym('x2');    % Underflow oil fraction, BetaU
x = [x1;x2];        
nx = 2;     % Number of states

%% System inputs

u = SX.sym('D');    % Input
du = SX.sym('du');      % Input change
s = SX.sym('s');        % Slack

%% System constraints

du0 = 0;

s_min = 0;
s_max = inf;
s0 = 0;

%% Model Equations

par = initHC();
Vo = par.VRF();     %
Vf = par.Vol_HC-par.VRF;
RhoO = 810;     % kg/m3

%% Model variables

x1 = SX.sym('x1');
x2 = SX.sym('x2');
x = [x1;x2];

%% Model inputs

u = SX.sym('u'); % Overflow valve opening

%% Experimental polynomial for Qsep

p2 = -5.5455 ;% (-1.322e+08, -5.671e+07)
p1 = 2.205;%  (6748, 1.13e+04)
p0 = 0.7209;%  (0.7378, 0.7919)

Qo = valveConstant*pi*0.0127^2*u*sqrt((2/RhoO)*(P2-par.Patm));

% Oil inflow
Qino = Betain.*Qin/2;

%% Polynomial is valid for [m3/h], therefore converting the units

QoPerHour = Qo*3600;
QinoPerHour = Qino*3600;

QsepPerHour = QinoPerHour.*(p2*QoPerHour.^2+p1*QoPerHour+p0);

Qsep = QsepPerHour/3600;

%% Model equations

xdot1 = (1/Vo)*(Qsep-x1*Qo); 
xdot2 = (1/Vf)*(Qino-Qsep-x2*Qu);

xdot = [xdot1;xdot2];

% Objective term
L = Q*s^2 + R*du^2;

Uu=[u;du;s];

ode = struct('x',x,'p',Uu,'ode',xdot,'quad',L);

%% Integrator

opts = struct('tf',h);
F = integrator('F', 'cvodes', ode, opts);

% preparing symbolic variable

w={};

preparing numeric variables and bounds
w0 = [];
lbw = [];
ubw = [];

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

% Formulate the NLP
 %initial state as variable
Xk = MX.sym('X0', nx);
w = {w{:}, Xk};
lbw = [lbw; x0(1); x0(2)];
ubw = [ubw; x0(1); x0(2)];
w0 = [w0; x0(1); x0(2)];

for k=0:ph-1
    % variable for the control in this control interval
    Uk = MX.sym(['U_' num2str(k)]);
    DUk = MX.sym(['DU_' num2str(k)]);
    Sk = MX.sym(['S_' num2str(k)]);
    Uk1 = [Uk; DUk; Sk];
    w = {w{:}, Uk1};
    
    lbw = [lbw; u_min; -du_max; s_min];
    ubw = [ubw; u_max; du_max; s_max];
    
    w0 = [w0; u0; du0; s0];
    
    
     if k <= mh
            if k == 0
                DUk = Uk - u0;
            else
                DUk = Uk(1) - Uk_1;
            end
            
            % Add input movement constraint
            g = {g{:}, DUk};
            lbg = [lbg; -du_max];
            ubg = [ubg; du_max];
            
      else 
            Uk(1) = Uk_1;
            DUk = 0;
      end
    
    % Integrate till the end of the interval
    
    Fk = F('x0',Xk,'p',Uk1);
    
    J = J + Fk.qf; 
    
    %loop
    Pk = Fk.xf;
    
    % States variables for state at end of interval
    % create states
    
    
    Xk = MX.sym(['X_' num2str(k+1)], nx);
    
    w = {w{:}, Xk};
    lbw = [lbw; 0; 0];
    ubw = [ubw; 1; 1];
    w0 = [w0; x0(1); x0(2)];
    
 %calculating delta U
       
   
    if k ~= ph - 1
            % for computing du
            Uk_1 = Uk1(1);
    end
    % Add continuinity constraint (closing shooting gap)
    g = {g{:}, Pk-Xk};
    lbg = [lbg; zeros(nx,1)];
    ubg = [ubg; zeros(nx,1)];
    
    % Soft constraint
    g = {g{:}, (Xk(2)-Sk-x2SP)};
    lbg = [lbg; -deltaX];
    ubg = [ubg; 0]; 
        
end

%% Create an NLP solver

prob = struct('f', J, 'x', vertcat(w{:}), 'g', vertcat(g{:}));

solver = nlpsol('solver', 'ipopt', prob);

%% Solve the NLP

sol = solver('x0', w0, 'lbx', lbw, 'ubx', ubw,'lbg', lbg, 'ubg', ubg);
w_opt = full(sol.x);


%% Various outputs

uOpt = w_opt(nx + 1);

obj = full(sol.f);

x1Horizon = w_opt(1:5:end);

P = w_opt(2:5:end);    % x2

M = w_opt(3:5:end);    % Control inputs

stats = solver.stats();    % information

iterations = stats.iter_count;  % Number of iterations performed

success = stats.success;    % 1 = success 0 = fail

exitFlag = stats.return_status();

QoValve = full(Qo);

firstx2 = P(2);

lastx2 = P(11);

timeElapsed = toc;  % ending timer

end

