Plot for concentration profile?

I have tried making code for RED, but am not sure about concentration profile and boundary conditions

%% 
clc;
clear;
close all;

% Given Inputs (Boundary Conditions)
Nm = 50;          % Number of membranes
Ns = 5;           % Number of stacks
Np = 20;          % Number of branches
Am = 0.06;        % Membrane Area (m^2)
T = 298.15;       % Temperature (K)
Vh = 11.989;      % Volume of high salinity stream (m^3)
Vl = 13.838;      % Volume of low salinity stream (m^3)
Ch0 = 61.35;      % Concentration of high salinity stream (mol/m^3)
Cl0 = 1.083;      % Concentration of low salinity stream (mol/m^3)
fRED = 0.5;       % Operability factor
Ra = 0.0011;      % Resistance anion membrane (ohm·m²)
Rc = 0.0011;      % Resistance cation membrane (ohm·m²)
Ll = 0.01;        % Thickness low salinity channel (m)
Lh = 0.02;        % Thickness high salinity channel (m)
a = 0.98;         % Selectivity
Sh = 0.0147;      % Conductivity high salinity (S/m)
Sl = 0.016;       % Conductivity low salinity (S/m)

% Constants
F = 96485;        % Faraday's constant (C/mol)
Rg = 8.314;       % Universal gas constant (J/mol·K)

% Calculations
B = Vl / Vh;

Ch = Ch0 - (B / (1 + B)) * fRED * (Ch0 - Cl0);
Cl = Cl0 + (1 / (1 + B)) * fRED * (Ch0 - Cl0);

Ve = 2 * Nm * Ns * (a * Rg * T / F) * log10(Ch / Cl);

q = a * (Vh / (Nm * Ns)) * F * (B / 1 + B) * fRED * (Ch0 - Cl0);

Rs = (Nm * Ns) / (Am * Np) * (Ra + Rc + (Ll / (Sl * Cl)) + (Lh / (Sh * Ch)));

ERED = (Rs / (Rs + 0.1 * Rs)) * Ve * q;
ERED_prime = ERED * 0.0000001 * 0.278;

% Display the Results
fprintf('Results:\n');

Results:

fprintf('Final High Salinity Concentration (Ch): %.4f mol/m^3\n', Ch);

Final High Salinity Concentration (Ch): 45.2046 mol/m^3

fprintf('Final Low Salinity Concentration (Cl): %.4f mol/m^3\n', Cl);

Final Low Salinity Concentration (Cl): 15.0711 mol/m^3

fprintf('Electromotive Force (Ve): %.4f V\n', Ve);

Electromotive Force (Ve): 6.0053 V

fprintf('Charge Transported (q): %.4f C\n', q);

Charge Transported (q): 315426.9292 C

fprintf('System Resistance (Rs): %.4f Ohm\n', Rs);

System Resistance (Rs): 15.3682 Ohm

fprintf('Power Output (ERED): %.4f W\n', ERED);

Power Output (ERED): 1722027.0339 W

fprintf('Power Output (ERED): %.4f W\n', ERED');

Power Output (ERED): 1722027.0339 W

% Performance Evaluation: ERED vs Membrane Area
Am_values = linspace(0.5*Am, 1.5*Am, 50);
ERED_values = zeros(size(Am_values));

for i = 1:length(Am_values)
    Am_var = Am_values(i);
    Rs_var = (Nm * Ns) / (Am_var * Np) * (Ra + Rc + (Ll / (Sl * Cl)) + (Lh / (Sh * Ch)));
    ERED_var = (Rs_var / (Rs_var + 0.1 * Rs_var)) * Ve * q;
    ERED_values(i) = ERED_var;
end

% Plotting
figure;
plot(Am_values, ERED_values, '-o', 'Color', 'b', 'MarkerSize', 5, 'MarkerFaceColor', 'b');
grid on;
xlabel('Membrane Area (m^2)');
ylabel('Power Output ERED (W)');
title('Performance Evaluation: ERED vs Membrane Area');

●

[Torsten]

There is no variation in the variable "ERED_var":

ERED_var = (Rs_var / (Rs_var + 0.1 * Rs_var)) * Ve * q = 1/1.1 * Ve * q,

and this value does not depend on Am.