function dhdt = TanksODE(x,u,d,p)
%% Parameters
N   = p(end);                                             % number of tanks
A   = p(1:N);                                                   % Tank area
Cv  = p(N+1:2*N+1);                                     % valve coefficient
DP  = p(2*N+2:3*N+2);                   % pressure drop over valve=constant
rho = p(end-1);                                                   % density

%% Variables
h = x;                                                         % tank level
z = u;                                           % input. valaves positions
%% ODE
q   = Cv.*z.*sqrt(DP./rho);                                         % flows

 for j=2:N+1
      if h(j-1)<1e-8
           q(j) = q(j-1);
      end
 end

dhdt = zeros(N,1);                                      % allocate memory
for i = 1 :N                                        % N+1= number of valves
     dhdt(i) = 1/A(i)*(q(i)-q(i+1));                    % tank mass balance     
end


