function dxdt = modelv2(x,u,d,p)
ME = x(1);    %[mol] mass of organic (extract) layer (acid, solvent, water)
MR = x(2);   %[mol] mass water (raffinate) layer (assumed to be pure water)
xAE = x(3);      % [-] mole fractions of acid in organic layer  (extract E)
xWE = x(4);     % [-] mole fractions of water in organic layer  (extract E)

F = u;                                        % [mol/s] Water and acid feed

S = d(1);                                    % [mol/s] organic solvent feed
xAF = d(2);                        % [-] mole frocation of acid in feed (F)

K = p(1);                  % [-] water acid equilibrium constant in extract
Kc1 = p(2);               % [1/s] proportional gain for holdup p-controller
Kc2 = p(3);     % [1/s] proportional gain for raffinate holdup p-controller

%% mathematical model 

E = Kc1*(ME+MR);
R = Kc2*MR;

dMEdt = (1+K)*xAF*F + S - E;

dMRdt = (1-(1+K)*xAF)*F - R;

dxAEdt = (1/ME)*( xAF(1-xAE*(1+K))*F - xAE*S);

dxWEdt = K*dxAEdt;

dxdt = [dMEdt;dMRdt;dxAEdt;dxWEdt];

end

