  function xend = Plant_three_tanks(x0,U,dt)

% calling CasADi
addpath('C:\Users\Alexander\Dropbox\Alexander\NTNU\Prosess-systemteknikk\MPC\CasADi')
import casadi.*

%% Parameter values
rho = 1000; % kg/m^3
g = 9.81; % m/s^2
A1 = 5; % m^2
A2 = 5; % m^2
A3 = 5; % m^2
Cv_0 = 1.7889*10^(-3);
Cv_1 = 3.2308*10^(-3);
Cv_2 = 2.1067*10^(-3);
Cv_3 = 3.2180*10^(-3);
DeltaP_0 = 50000; % [Pa], constant pressure drop over z0

% Nominal levels
h10 = 1.0;
h20 = 1.2;
h30 = 0.7;   
%% Model variables
h1 = SX.sym('h1'); 
h2 = SX.sym('h2');
h3 = SX.sym('h3');
x = [h1;h2;h3];

%% System inputs
z0 = SX.sym('z0');
z1 = SX.sym('z1');
z2 = SX.sym('z2');
z3 = SX.sym('z3');
u = [z0;z1;z2;z3];

%%  Dynamic equations
    % Nonlinear
%     xdot = [1/A1*( Cv_0*U(1)*sqrt(DeltaP_0)   - Cv_1*U(2)*sqrt(rho*g*x(1)) );
%             1/A2*( Cv_1*U(2)*sqrt(rho*g*x(1)) - Cv_2*U(3)*sqrt(rho*g*x(2)) );
%             1/A3*( Cv_2*U(3)*sqrt(rho*g*x(2)) - Cv_3*U(4)*sqrt(rho*g*x(3)) )];
    % Linear
    xdot = [1/A1*( Cv_0*U(1)*sqrt(DeltaP_0)  - Cv_1*U(2)*sqrt(rho*g*h10) );
            1/A2*( Cv_1*U(2)*sqrt(rho*g*h10) - Cv_2*U(3)*sqrt(rho*g*h20) );
            1/A3*( Cv_2*U(3)*sqrt(rho*g*h20) - Cv_3*U(4)*sqrt(rho*g*h30) )];
%% Integrating the system
% Formulate discrete time dynamics
% CVODES from the SUNDIALS suite
ode = struct('x',x,'p',u,'ode',xdot);

% Building the integrator
opts = struct('tf',dt);
F = integrator('F', 'cvodes', ode, opts);

sim = F('x0',x0,'p',U);

xend = full(sim.xf);
