function dxdt=uebung_3_system(t,x,p)

%Zerlegung des Vektors
c_A = x (0*p.n+1 : 1*p.n);
c_B = x (1*p.n+1 : 2*p.n);
theta_A = x (2*p.n+1 : 3*p.n);
theta_B = x (3*p.n+1 : 4*p.n);

%Randbed.

if t >= 0 & t <= 0.1
    c_A_in = 1e-5;  % [mol / m³]
    c_B_in = 1e-5;  % [mol / m³]
    
else  
    c_A_in = 0;
    c_B_in = 0;
end   
% Hilfsvariable

r_adde_A(:,1) = p.k_ad_A*(c_A .* (1- (theta_A + theta_B)) - 1./p.k_adde_A*theta_A);
r_adde_B(:,1) = p.k_ad_B*(c_B .* (1- (theta_A + theta_B)) - 1./p.k_adde_B*theta_B);

% Ableitung der Konzentration

dc_Adt(1,1) = p.D_A * (( c_A(2) - c_A(1))/p.dz^2 - (c_A(1) - c_A_in)/p.dz^2*2) ...
    -p.v*( c_A(1) - c_A_in ) /p.dz -r_adde_A(1)*p.a;  % 1.Kontrollvol

dc_Adt(p.n,1) = p.D_A * ( 0 - (c_A(p.n) - c_A(p.n-1))/p.dz^2 ) ...
    -p.v*(c_A(p.n) - c_A(p.n-1))/p.dz - r_adde_A(p.n)*p.a; % letztes Kontrollvol

dc_Adt(2:p.n-1,1) =  p.D_A * ( (c_A(3:p.n) - c_A(2:p.n-1))/p.dz^2 - (c_A(2:p.n-1)-c_A(1:p.n-2))/p.dz^2)...
    -p.v*(c_A(2:p.n-1) - c_A(1:p.n-2))/p.dz - r_adde_A(2:p.n-1)*p.a; % mittl. Kontrollvol


% für B

dc_Bdt(1,1) = p.D_B * (( c_B(2) - c_B(1))/p.dz^2 - (c_B(1) - c_B_in)/p.dz^2*2) ...
    -p.v*( c_B(1) - c_B_in ) /p.dz -r_adde_B(1)*p.a;  % 1.Kontrollvol

dc_Bdt(p.n,1) = p.D_B * ( 0 - (c_B(p.n) - c_B(p.n-1))/p.dz^2 ) ...
    -p.v*(c_B(p.n) - c_B(p.n-1))/p.dz - r_adde_B(p.n)*p.a; % letztes Kontrollvol

dc_Bdt(2:p.n-1,1) =  p.D_B * ( (c_B(3:p.n) - c_B(2:p.n-1))/p.dz^2 - (c_B(2:p.n-1)-c_B(1:p.n-2))/p.dz^2)...
    -p.v*(c_B(2:p.n-1) - c_B(1:p.n-2))/p.dz - r_adde_B(2:p.n-1)*p.a; % mittl. Kontrollvol


dtheta_Adt = r_adde_A/p.c;
dtheta_Bdt = r_adde_B/p.c;

dxdt = [dc_Adt; dc_Bdt; dtheta_Adt; dtheta_Bdt];


