How to calculate initial concentration of first-order reaction in CSTR?

(This question comes from a Chemical Engineering background, but basic Math and Programming knowledge is also required. I hope it still falls in the scope of engineering.stackexchange.com If not, please move.)

Let's assume a (slow hydrolysis) first-order reaction of a component D in a CSTR, which is described by the well-known ODE:

$$\frac{d[D]}{dt} = \frac{\dot{V}_{in}}{V_{liq}}([D]_{in} - [D]) - k_1 [D]$$

where $$[D]$$ is concentration of D, $$\dot{V}_{in} = q_{in}$$ is input volumetric flow, $$V_{liq}$$ is the reactor volume, and $$k_1$$ is the first order rate constant. The reactant is dissolved in water and converted to two gaseous components whose solubility is practically 0. Thus, the formation rate is proportional to the reaction rate of D, $$k_1[D]$$.

I have put this in a small Matlab/Octave (main) script:

% Define basic parameters
V_liq   = 2500;     % m^3
k_1     = 0.320;    % 1/d
D_in    = 180;      % kg/m^3
q_in    = 40;       % m^3/d
Y_M     = 420;      % m^3/t [D]
Y_N     = 389;      % m^3/t [D]
D_steady    = (D_in*q_in/V_liq)/(q_in/V_liq + k_1);       % kg/m^3

[D_steady, q_gas_M/24, q_gas_N/24] % Second element should = 120 (m^3/h)

ans =

8.5714 120.0000 111.1429

Now, I am trying to validate my dynamic simulation by passing a constant feed vector. The feed occurs every hour, with the fresh liquid replacing the same amount of reaction liquid, thus maintaining $$V_{liq}$$ constant. (Strictly speaking, this is not a CSTR anymore, rather a semi- or quasi-CSTR.)

I have put the model equations in separate function (using a for-loop for the feed):

function resultsArray = reactorModel1(D_feed, initialValue, parameters, D_in)
% Simulate gas production per hour with algebraic reactor model
% Feed is solved via a for-loop

hoursTotal  = length(D_feed);
k_1         = parameters(1);
Y_M         = parameters(2);
Y_N         = parameters(3);
V_liq       = parameters(4);
resultsArray = zeros(hoursTotal, 3);
t           = 1/24;

liquid_feed = D_feed/(D_in*1e-3);

initialValue4Model0 = (initialValue*(V_liq - liquid_feed(1))*1e-3 ...
+ D_feed(1))*1e3/V_liq; % kg/m^3
resultsArray(1, 1) = initialValue4Model0*exp(-k_1*t);
% Simple for-loop with feed as vector per hour
for pHour = 2:hoursTotal
initialValue4Model = (resultsArray(pHour-1, 1)*(V_liq - liquid_feed(pHour))*1e-3 ...
+ D_feed(pHour))*1e3/V_liq; % kg/m^3
resultsArray(pHour, 1) = initialValue4Model*exp(-k_1*t);
end
resultsArray(:, 2) = V_liq*Y_M*k_1*resultsArray(:, 1)*1e-3; % m^3/d
resultsArray(:, 3) = V_liq*Y_N*k_1*resultsArray(:, 1)*1e-3; % m^3/d
end

Which is called in the second part of my main script, where I (wrongfully?) assume $$[D](t=0) = [D]_{steady}$$:

%% Part 2: Dynamic calculation with constant feeding
parameters  = [k_1, Y_M, Y_N, V_liq];
t_max       = 168; % hours
D_feed      = q_in*D_in/24*1e-3*ones(t_max, 1);
resultsArray = reactorModel1(D_feed, initialValue, parameters, D_in);
[resultsArray(end, 1), resultsArray(end, 2:3)/24]

ans =

8.5223 119.3119 110.5056

As you can see, all three components at $$t_{max}$$ are lower than in steady-state. If you take a look at the progression of D (figure; plot(resultsArray(:, 1))), you see it slightly decreasing over time. However, all parameters (liquid volume, kinetic constant, flowrate of feed, concentration of D in water) are the same between steady-state and dynamic simulation!

So my question is: Is there a way to **analytically determine $$[D](t=0)$$ ** so that $$[D](t)$$ becomes constant? (As stated above in the comment, I need q_gas_M = 120, since only gaseous flows are continuously measured.)

(Partly answering my own question) It appears that if I use a different solution method, I actually get the expected result for constant input.

The alternative solution consists in solving the whole ODE (including mass transport terms) symbolically, and then modifying the resulting algebraic equation to accept t and q_in as vectors, rather than just scalars.

function resultsArray = reactorModel2(D_feed, initialValue, parameters, D_in)
% Simulate gas production per hour with algebraic reactor model
% Feed is solved as part of an algebraic equation (symbolic solution of CSTR ODE)

hoursTotal  = length(D_feed);
t           = (1:hoursTotal).'/24;
q_in        = D_feed/D_in*1e3*24; % m^3/h

k_1         = parameters(1);
Y_M         = parameters(2);
Y_N         = parameters(3);
V_liq       = parameters(4);
resultsArray = zeros(hoursTotal, 3);
% Symbolic solution of CSTR ODE
resultsArray(:, 1) = (D_in*q_in + exp(-(t.*(q_in + V_liq*k_1))/V_liq).*(q_in*(initialValue - D_in) ...
+ V_liq*k_1*initialValue))./(q_in + V_liq*k_1); % kg/m^3

resultsArray(:, 2) = V_liq*Y_M*k_1*resultsArray(:, 1)*1e-3; % m^3/d
resultsArray(:, 3) = V_liq*Y_N*k_1*resultsArray(:, 1)*1e-3; % m^3/d
end

The corresponding call from the console

resultsArray2 = reactorModel2(D_feed, initialValue, parameters, D_in);
[resultsArray2(end, 1), resultsArray2(end, 2:3)/24]

gives

ans =

8.5714 120.0000 111.1429

When I compare both solutions graphically, you can see that they slightly differ. Not sure why. However, as soon as I start using a dynamic input vector (maintaining the overall sum of fed amount constant), Method 1 gives far more plausible results than Method 2. But this should be adressed in another question, I guess.