For loop gives the same output

I have a problem with this code that it gives the same output when the loop reaches the 4th column of the inputs, the loop gives the same results (when it should not). For the hell of me I don't know what is wrong! It has to be in the loop not in the numbers! I ran them indvidually for all scenarios and got unique output for every scenario!

%% Defining Scenario Independent Paramters:
%For Rocket
%Payload mass:
M_L=1;%kg
%Number of fins:
N=3;
%Shell density:
rho_s=2700; %kg/m^3
%Propellant density:
rho_p=1772; %kg/m^3
%Shell working stress:
sig_s=60*10^6; %pa
%Gravity:
g=9.81; % m/s^2

%For air:
%Atmorspheric TeM_perature at sea level:
T_atm=298;
%Specific heat ratio of air:
gamma=1.4;
%Density of air at sea level:
rho=1.225;
%Pressure at sea level:
P_a=101325; %pa

%%My 9 scenarios that are defined the text data are:
%Maximum Altitude, h_max= 0.3048*[20000 20000 20000 2000 2000 10000 30000 10000 30000]; %m
%Normalized max acceleration, a_max= [10 10 10 5 20 10 10 5 20]
% Static margin , SM = [1 2 3 2 2 2 2 1 3]
%% Calculating D_max and L_max through the 9 secenarios
Data=readmatrix('Data');% Each col is 1 test
D_max=zeros(1,9);
L_max=zeros(1,9);
for loop=1:9
    h_max=Data(1,loop);
    a_max=Data(2,loop);
    SM=Data(3,loop);
    R_max=1+a_max;
    W_eq=sqrt((h_max*g)/(((log(R_max)/2)*(log(R_max)-2))+((R_max-1)/R_max)));
    t_bmax=(R_max-1)*W_eq/(g*R_max);
    M_eq=W_eq/sqrt(gamma*287*T_atm);
    P_c=P_a*(1+(((gamma-1)/2)*M_eq^2))^(gamma/(gamma-1));
    P_0_a=P_c/P_a;
    
    %Iteration for best L, D
    % Creating Length and Diameter Limiatation for Iteration
    L=linspace(0.0001,4,5000); %m
    D=linspace(0.0001,2,5000); %m
    
    %Iteration through all the parameters:
    for i=1:length(D);
        for j=1:length(L)
            delta= (P_c/(2*sig_s))*D(i); % thickness of shell m
            M_n=delta*rho_s*pi*D(i)*(D(i)+sqrt(D(i)^2+(D(i)^2/4)));
            M_f=((D(i)^2)/2)*delta*rho_s;
            M_f_b=(pi*D(i)*rho_s*D(i)*delta);
            M_s=(pi*D(i)*rho_s*L(j)*delta)+M_n+M_f+M_f_b;%%%%%%%%%%
            M_p=(R_max-1)*(M_s+M_L);
            L_p=(M_p/(pi*D(i)^2*rho_p/4));
            if L_p<L(j)+D(i)
                % Center of Pressure for Rocket Nose
                X_n=(2/3)*D(i);%m
                CN_n=2;
                
                %Center of Pressure for Rocket Fin
                a=D(i);%m
                s=D(i);%m
                b=0;
                m=a-b;
                fin_hyp=sqrt(2)*D(i);%m
                X_f=D(i)+L(j);%m
                delta_X_f=((m*(a+2*b))/(3*(a+b)))+((1/6)*(a+b-((a*b)/(a+b))));%m
                X_f_dash=X_f+delta_X_f;%m
                CN_f=(4*N*(s/D(i))^2)/(1+sqrt(1+((2*fin_hyp)/(a+b))^2));%m
                k_fb=1+((D(i)/2)/(s+(D(i)/2)));%m
                CN_fb=k_fb*CN_f; %m;
                X_cp=((CN_n*X_n)+(CN_fb*X_f_dash))/(CN_n+CN_fb);
                
                
                %Center of Gravity for Rocket Cone:
                X_ncg=2*D(i)/3;
                M_n=delta*rho_s*pi*D(i)*(D(i)+sqrt(D(i)^2+(D(i)^2/4)));
                M_c=M_L+M_n;
                
                %Center of Gravity for Rocket Fin
                X_f_cg=(2*D(i)/3)+L(j)+D(i);
                M_f=((D(i)^2)/2)*delta*rho_s;
                
                %Center of Gravity for Rocket Tube 1
                L_1=L(j)+D(i)-L_p;
                xL_1_cg=(L_1/2)+D(i);
                M_L_1=pi*D(i)*L_1*delta*rho_s;
                
                %Center of Gravity for Rocket Tube 2
                L_2=L_p;
                xL_2_cg=(L_2/2)+D(i)+L_1;
                M_L_2=(pi*D(i)*L_2*delta*rho_s)+M_p;
                X_cg=((M_c*X_ncg)+(M_L_1*xL_1_cg)+(xL_2_cg*M_L_2)+(X_f_cg*M_f*3))/(M_c+M_L_1+M_L_2+(M_f*3));
                
                cond(i,j)=X_cp-X_cg-(D(i)*SM);
                if cond(i,j)<0 || cond(i,j)>0.00001
                    continue
                end
                L_D(i,j)=L(j)/D(i);
                lamda(i,j)=M_L/(M_s+M_p);
            else
                continue
            end
        end
        
    end
    
    %Output
    % Matrix indexing for the D_max and L_max
    [M,I] = max(lamda,[],'all','linear');
    [row,col] = ind2sub(size(lamda),I);
    D_max(1,loop)=D(row)
    L_max(1,loop)=L(col)
end

Answer 1

[Bjarke Skogstad Larsen]

You need to initialize the matrixes: L_D, lamda and cond.

Always initialize variables before a loop when you only set part of the variable at a time inside the loop, like

lamda(i,j)=M_L/(M_s+M_p);

needs to be initialized like so:

    lamda = nan(length(D),length(L));
    for ii=1:length(D)

I will include the full working code below (I changed i to ii and j to jj as I don't like using these variables in Matlab due to the possibility of complex values):

clc
clear
%% Defining Scenario Independent Paramters:
%For Rocket
%Payload mass:
M_L=1;%kg
%Number of fins:
N=3;
%Shell density:
rho_s=2700; %kg/m^3
%Propellant density:
rho_p=1772; %kg/m^3
%Shell working stress:
sig_s=60*10^6; %pa
%Gravity:
g=9.81; % m/s^2

%For air:
%Atmorspheric TeM_perature at sea level:
T_atm=298;
%Specific heat ratio of air:
gamma=1.4;
%Density of air at sea level:
rho=1.225;
%Pressure at sea level:
P_a=101325; %pa

%%My 9 scenarios that are defined the text data are:
%Maximum Altitude, h_max= 0.3048*[20000 20000 20000 2000 2000 10000 30000 10000 30000]; %m
%Normalized max acceleration, a_max= [10 10 10 5 20 10 10 5 20]
% Static margin , SM = [1 2 3 2 2 2 2 1 3]
%% Calculating D_max and L_max through the 9 secenarios
Data=readmatrix('Data');% Each col is 1 test
D_max=zeros(1,9);
L_max=zeros(1,9);
for loop=1:9
    disp(num2str(loop));
    h_max=Data(1,loop);
    a_max=Data(2,loop);
    SM=Data(3,loop);
    R_max=1+a_max;
    W_eq=sqrt((h_max*g)/(((log(R_max)/2)*(log(R_max)-2))+((R_max-1)/R_max)));
    t_bmax=(R_max-1)*W_eq/(g*R_max);
    M_eq=W_eq/sqrt(gamma*287*T_atm);
    P_c=P_a*(1+(((gamma-1)/2)*M_eq^2))^(gamma/(gamma-1));
    P_0_a=P_c/P_a;
    
    %Iteration for best L, D
    % Creating Length and Diameter Limiatation for Iteration
    L=linspace(0.0001,4,5000); %m
    D=linspace(0.0001,2,5000); %m
    
    %Iteration through all the parameters:
    L_D = nan(length(D),length(L));
    lamda = nan(length(D),length(L));
    cond = nan(length(D),length(L));
    for ii=1:length(D)
        for jj=1:length(L)
            delta= (P_c/(2*sig_s))*D(ii); % thickness of shell m
            M_n=delta*rho_s*pi*D(ii)*(D(ii)+sqrt(D(ii)^2+(D(ii)^2/4)));
            M_f=((D(ii)^2)/2)*delta*rho_s;
            M_f_b=(pi*D(ii)*rho_s*D(ii)*delta);
            M_s=(pi*D(ii)*rho_s*L(jj)*delta)+M_n+M_f+M_f_b;%%%%%%%%%%
            M_p=(R_max-1)*(M_s+M_L);
            L_p=(M_p/(pi*D(ii)^2*rho_p/4));
            if L_p<L(jj)+D(ii)
                % Center of Pressure for Rocket Nose
                X_n=(2/3)*D(ii);%m
                CN_n=2;
                
                %Center of Pressure for Rocket Fin
                a=D(ii);%m
                s=D(ii);%m
                b=0;
                m=a-b;
                fin_hyp=sqrt(2)*D(ii);%m
                X_f=D(ii)+L(jj);%m
                delta_X_f=((m*(a+2*b))/(3*(a+b)))+((1/6)*(a+b-((a*b)/(a+b))));%m
                X_f_dash=X_f+delta_X_f;%m
                CN_f=(4*N*(s/D(ii))^2)/(1+sqrt(1+((2*fin_hyp)/(a+b))^2));%m
                k_fb=1+((D(ii)/2)/(s+(D(ii)/2)));%m
                CN_fb=k_fb*CN_f; %m;
                X_cp=((CN_n*X_n)+(CN_fb*X_f_dash))/(CN_n+CN_fb);
                
                
                %Center of Gravity for Rocket Cone:
                X_ncg=2*D(ii)/3;
                M_n=delta*rho_s*pi*D(ii)*(D(ii)+sqrt(D(ii)^2+(D(ii)^2/4)));
                M_c=M_L+M_n;
                
                %Center of Gravity for Rocket Fin
                X_f_cg=(2*D(ii)/3)+L(jj)+D(ii);
                M_f=((D(ii)^2)/2)*delta*rho_s;
                
                %Center of Gravity for Rocket Tube 1
                L_1=L(jj)+D(ii)-L_p;
                xL_1_cg=(L_1/2)+D(ii);
                M_L_1=pi*D(ii)*L_1*delta*rho_s;
                
                %Center of Gravity for Rocket Tube 2
                L_2=L_p;
                xL_2_cg=(L_2/2)+D(ii)+L_1;
                M_L_2=(pi*D(ii)*L_2*delta*rho_s)+M_p;
                X_cg=((M_c*X_ncg)+(M_L_1*xL_1_cg)+(xL_2_cg*M_L_2)+(X_f_cg*M_f*3))/(M_c+M_L_1+M_L_2+(M_f*3));
                
                cond(ii,jj)=X_cp-X_cg-(D(ii)*SM);
                if cond(ii,jj)<0 || cond(ii,jj)>0.00001
                    continue
                end
                L_D(ii,jj)=L(jj)/D(ii);
                lamda(ii,jj)=M_L/(M_s+M_p);
            else
                continue
            end
        end
        
    end
    
    %Output
    % Matrix indexing for the D_max and L_max
    [M,I] = max(lamda,[],'all','linear');
    [row,col] = ind2sub(size(lamda),I);
    D_max(1,loop)=D(row)
    L_max(1,loop)=L(col)
end

[Mohsen Nashel]

Thank you so much! It worked!