Unknown For Loop Placement Error / Increment Error

Hello, I have a script for modeling a CO2 Scrubbing system found in smoke stacks of power plants which computes the concentration of CO2 using the finite difference method. The script iteratively computes the concentration at a (z,r) coordinate and compares it to a local index of (i,j). I am being asked to compute the concentration of the exhaust gas for a range of smoke stack heights. The script successfully solves for the concentration when individual values of length are used, however, I would like to solve for 50 different concentrations with 50 different lengths. To do this, I took the section of code which workes for individual values, and encased it in a for loop which updates the length after the concentration gradient has been determined. My problem is when running the script, it appears that only the first value of the length is used, yet it still returns 50 values of concentrations (they are all the same since the first length is being utilized). I would like to know if anyone can recognize where the issuse might occur, potentially from loop placement / placement of where the length value is incrimented.

***Additionally, the script is specifically looking to calculate the concentration at the exit of the smoke stack which is represented as C_TOP. This does not cause addtional problems with computation, but can be used with C_avg to see that the results are only evaluate for L(ii)=L(1).

Thanks for any help,

Kyle

%% System inputs %%

%Operational Parameters:
C0 = 0.3;        %Inlet mass fraction [-]
Cw = 0.0;        %Reference mass fraction with amine [-]
kw = 1;          %Convective mass transfer coefficent [m/s]
Q  = 1;          %Volumetric flow rate [m3/s]
D  = 0.02;       %Diffusion coefficent (CO2 - amine) [m2/s]


%Design Parameters:
L = linspace(0.1,20,50);        %length [m] between 0 and 20 meters
R = 0.8;         %radius [m] between 0.8 and 1.6 meters
ax = Q/(pi*R^2); %Representative advective velocity

%Design Assessment Parameters:
C_target = 0.1;     %target concentration
Vs0      = 20;      %cost of scrubber + operational costs
ms       = 1.4;
mc       = 100;



%% Model Solution %% 

% numerical:
N = 50;        %number of nodes
Nr = N;        %number of nodes along r
Nz = N;        %number of nodes along z

% tolerance tolerance for solution
tolerance = 1.0e-4;
% define a 2D grid with nodal values of variable:
C = zeros( Nz, Nr, N );

% set up: ---------------------------------------------------------------
for ii = 1: N       %start of for loop used to incrment L
    
%everything from here down works for individual values of L and R but in
%that case L = L, not L = L(ii)

% spacing between nodes in each direction:
deltar = R / ( Nr - 1 );
deltaz = L(ii) / ( Nz - 1 );       %this is where the length value is update for each itteration of for loop (or so i think)

% solution procedure: ---------------------------------------------------
Cold(:,:,ii)    = C(:,:,ii);       % auxiliary array to store solution
error   = 1.0e+8;                  % initialize error as a large number
nite    = 0;                       % iteration counter
nitemax = 1e+5;                    % maximum number of iterations

    % iterate until convergence:
    while ( error > tolerance ) && ( nite < nitemax )

    % update iteration counter:
    nite = nite + 1;

    % estimate new value in nodes from the old values:
        for i = 1 : Nz
            for j = 1 : Nr
            
            % define coordinates for point with indexes (i, j):
            r = deltar * ( j - 1 );
            z = deltaz * ( i - 1 );
            
            % update material properties:
            % - if needed, k, S can be updated as function of r, z, Told
 
            % handle boundaries:
                if     i ==  1 % bottom
                C( i, j, ii ) = C0;
                        
                elseif i == Nz % top
                C( i, j , ii ) =  Cold( i-1, j ); 
                        
                elseif j ==  1 % left 
                C( i, j, ii ) =  Cold( i, j+1 );
                            
                        
                elseif j == Nr % right
                C( i, j, ii ) = ( Cold( i, j-1, ii ) + ( kw  * deltar / D ) * Cw  ) / ...
                            ( 1              + ( kw  * deltar / D ) );

                else % we are inside the domain:

                % coefficients:               
                Cleft   = -(D/(2*r*deltar)) - (D/(deltar^2));
                Cright  =  (D/(2*R*deltar)) - (D/(deltar^2));
                Cbottom = -(ax/(2*deltaz))  - (D/(deltaz^2));
                Ctop    =  (ax/(2*deltaz))  - (D/(deltaz^2));
                Ccenter =  Cbottom + Ctop +  Cleft  + Cright;
                
                % update solution:
                C( i, j, ii ) = ( Cleft   * Cold( i  , j-1, ii ) + ...
                              Cright  * Cold( i  , j+1, ii ) + ...
                              Cbottom * Cold( i-1, j, ii   ) + ...
                              Ctop    * Cold( i+1, j, ii   ) + ...
                              0 ) / Ccenter;
                          
                end
            end
        
        end
    
        % update error:
        error = norm( C(:,:,ii) - Cold(:,:,ii) );
    
        % display convergence:
        %fprintf('nite = %i, error = %2.2e \n ', error );

        % store the solution from this iteration:
        Cold(:,:,ii) = C(:,:,ii);

   end
end

%% Calculating Cavg & CO2 Absorbed %%

%extracting C values at top of pipe
C_TOP = C(N,:,:);

%Cavg calculation
syms r
C_avg = zeros(50,1);

for i = 1:length(C_TOP)

    C_avg(i) = ((1/(pi*R^2))*int(2*pi*r,r,0,R))*sum(C_TOP(:,:,i))/N;

end

C_avg