function dxdt = CL_neutral_power(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

% Compute neutral power assumed to be plant output in the current state 
Salpha = sources(abs(T),0,p);   % note that the alpha power doesn't depend on the auxiliary power, so I put 0
P_neutron=4*Salpha.S_alpha*V;

%% Linear controller
e = C_ss.P_neutron - P_neutron; % control error: setpoint is P_neutron_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