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
