function dTdt = Tank(x,u,p,d)
% model for two tanks in series with heating element in tank 1

%parameters
V1  = p(1);                                            % [m3] tank 1 volume
V2  = p(2);                                            % [m3] tank 2 volume
rho = p(3);                                               % [m3/kg] density
cp  = p(4);                                    % [kJ/kg/degC] specific heat


%states: 2 =>[T1;T2]
T = x;                                    % [degC] tank 1 and 2 temperature

% input
q = u;                                                        % [m3/s] flow

%  disturbances
T0 = d(1);                                       % [degC] inlet temperature
TD = d(2);               % [degC] inlet temperature of flow added to tank 2
qD = d(3);                              % [m3/s] flow rate added to  tank 2
Q  = d(4);                                                % [kW] heat input

% energy balance. ODE
dTdt(1) = 1/V1*(q*(T0-T(1)) + Q/(rho*cp));                         % tank 1
dTdt(2) = 1/V2*(q*(T(1)-T(2))+qD*(TD-T(2)));                       % tank 2

end