function ode = TESmodel(x,u,d,param)

%% Parameters
V   = param(1);                                               % tank volume
rho = param(2);                                                   % density
p0  = param(3);                                        % reference pressure
m0  = param(4);                                    % reference network mass
eps = param(5);                                     % compressibility coeff

%% States
Vh  = x(1);                                              % hot layer volume
p   = x(2);                                              % network pressure

%% Disturbances
qS   = d(1);                                      % hot waer flowrate demand
q    = d(2);
%% Inputs
qSP = u(1);                                 % hot water direct to consumers
qTP = u(2);
qEB = u(3);                                % hot water from electric boiler

if   Vh  >= V
     qD  = qS-qSP; % dumped
     qST = 0;      % discharge
elseif Vh< V
     qD  = 0;       % dump
     qST = qS-qSP;  % discharge
end



%% Hot layer tank balance 
dVhdt = qST-qTP;                                             % mass balance

%% Pipeline balance
dpdt = p0/m0/eps*rho*(qSP+qTP+qEB-q); 

%% Function output
ode = [dVhdt;dpdt];                                % Differential equations

