function dxdt = CL_thermal_energy_balance(t,x,time,C_ss,p)

q_e = 1.60217657E-19; % electron charge [coulomb] or Boltzmann constant [J/eV]
V = 2*pi^2*p.kappa*p.R_0*p.a^2; % plasma volume [m^3]
T = x(1); % plant state
z = x(2:end); % controller state

%% Linear controller
e = C_ss.T_0 - T; % control error: setpoint is T_0
dzdt =   C_ss.A*z + C_ss.B*e;
dP_aux = C_ss.C*z + C_ss.D*e;
S_aux = (dP_aux + C_ss.P_0)/V;

%% Nonlinear plant
S = sources(abs(T),dP_aux+C_ss.P_0,p);
dT_eVdt = 1/(3*p.n*q_e)*(S.S_alpha+S.S_ohm+S_aux-S.S_rad-S.S_cond);

%% Total system time derivative
dxdt = [dT_eVdt; dzdt];
end