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

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

%% Disturbances
q   = d(1);                                      % hot waer flowrate demand

%% Inputs
qST = u(1);                                       % hot water charging tank
qTP = u(2);                                         % hot water discharging
qSP = u(3);                                 % hot water direct to consumers
qEB = u(4);                                % hot water from electric boiler

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

%% 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

