function dxdt = TanksODE(x,u,d,p)
% energy balance for heated tanks in series

%% Parameters
F  = p(1);                                                       % flowrate
m1 = p(2);                                               % holdup in tank 1
m2 = p(3);                                               % holdup in tank 2
cp = p(4);                                                  % specific heat

%% States
T1 = x(1);                                          % temperature in tank 1
T2 = x(2);                                          % temperature in tank 2

%% Input
Q1 = u;                                                    % heat to tank 1

%% Disturbances
T0 = d(1);                                              % inlet temperature
Q2 = d(2);                                                 % heat to tank 2

%% Model. Energy balance. Assume constant volume
dT1dt = F/m1*(T0-T1) + Q1/(m1*cp);                                 % tank 1 
dT2dt = F/m2*(T1-T2) + Q2/(m2*cp);                                 % tank 2

dxdt  = [dT1dt; dT2dt];                                     % system of ODE