function  [output] = SysThreeTanks(h,u,p,F)
% Model parameters

rho = p(1); % kg/m^3    
A1= p(2); % m^2
A2= p(3); % m^2
A3= p(4); % m^2
g = 9.81; % m/s^2

% Disturbances
F4 = F(1);
F5 = F(2);
F6 = F(3);
F7 = F(4);
F8 = F(5);
F9 = F(6);

% Valve positions
z1=u(1);
z2=u(2);
z3=u(3);
z0=u(4);

% Valve constants
Cv_0 = 1.7889*10^(-3);
Cv_1 = 3.2308*10^(-3);
Cv_2 = 2.1067*10^(-3);
Cv_3 = 3.2180*10^(-3);

% With constant pressure drop over valves
q0 = Cv_0*z0*sqrt(50000);
q1 = Cv_1*z1*sqrt(rho*g*1);
q2 = Cv_2*z2*sqrt(rho*g*1.2);
q3 = Cv_3*z3*sqrt(rho*g*0.7);

% Differential equations
dxdt(1) = 1/A1*(q0-q1 + F4-F5);
dxdt(2) = 1/A2*(q1-q2 + F6-F7);      
dxdt(3) = 1/A3*(q2-q3 + F8-F9);
% q = [q0,q1,q2,q3];
% Making a vertical vector with states. 
output = [dxdt]';
end

