Please check with the following
clear, clc
V0 = 50; %initial volume (m^3)
c0 = 0; %initial concentration
Qin = 1.25 * 60; %inflow (m^3/hour)
Qout = 1.1 * 60; %outflow (m^3/hour)
cin = 100; %pollutant concentration (mg/L)
t_final = 50; %hours
dt = 0.5; %time step
T = [dt:0.5:50];
t = 0;
%
k = 0;
%
%Initialize Loop Variables
cold = c0;
vold = V0;
inflow=1; % inflow status variable
for n = [1 : (t_final/dt)]
t = t + dt;
k = 0;
if vold>=200 || inflow==0
inflow=0; % say inflow stopped
if vold<=150
inflow=1; % say inflow start
end
end
if inflow==0
Qin=0;
else
Qin=1.25 *60;
end
alpha = (Qin / (vold + (Qin - Qout) * t)) + k;
beta = (Qin*cin)/(vold+(Qin - Qout) * t);
c(n) = cold - dt * alpha * cold + dt * beta;
V(n) = vold + (Qin - Qout) * dt;
cold = c(n);
vold = V(n);
end
yyaxis left
plot(T,V)
ylabel('Volume (m^3)')
xlabel('Time (hours)')
yyaxis right
plot(T,c)
ylabel('Concentration (mg/L)')
disp(' t V c')
format compact
disp([T' V' c'])
more(10)
