diff --git a/Dipolar-Gas-Simulator/+Scripts/createPhaseDiagram.m b/Dipolar-Gas-Simulator/+Scripts/createPhaseDiagram.m
index 6097ca9..0f5e796 100644
--- a/Dipolar-Gas-Simulator/+Scripts/createPhaseDiagram.m
+++ b/Dipolar-Gas-Simulator/+Scripts/createPhaseDiagram.m
@@ -1,7 +1,7 @@
 function createPhaseDiagram()
     % Define axis values
-    SCATTERING_LENGTH_RANGE = [79.0 80.0 81.0 82.0 83.0 84.0 85.0 86.0 87.0 88.0 89.0 90.0];
-    NUM_ATOMS_LIST = 50000:5000:105000;
+    SCATTERING_LENGTH_RANGE = [80.00 81.04 82.08 83.12 84.17 85.21 86.25 87.29 88.33 89.38 90.42 91.46 92.50 93.54 94.58 95.62 96.67 97.71 98.75 99.79 100.83 101.88 102.92 103.96 105.00];
+    NUM_ATOMS_LIST = [100000 304167 508333 712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 4591667 4795833 5000000];
     
     xlen = length(NUM_ATOMS_LIST);
     ylen = length(SCATTERING_LENGTH_RANGE);
diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m
index b06b81b..ef14eef 100644
--- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m
+++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m
@@ -548,56 +548,40 @@ Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
 SaveDirectory = './Results/Data_3D/AnisotropicTrap/TiltedDipoles45';
 JobNumber     = 0;
 Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
-%%
-SaveDirectory = './Results/Data_3D/GradientDescent/aS_088_theta_000_phi_000_N_750000';
-JobNumber     = 0;
-Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
-%%
-SaveDirectory = './Results/Data_3D/PhaseDiagram/ImagTimePropagation/aS_9.458000e+01_theta_000_phi_000_N_304167';
-JobNumber     = 0;
-Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
-%%
-% Parameters you can set before the loop
-N         = 100000;
-theta     = 20.0;
-phi       = 0;
-JobNumber = 0;
-
-% Values of aS to loop over
-aS_values = [79.0 80.0 81.0 82.0 83.0 84.0 85.0 86.0 87.0 88.0 89.0 90.0];
-
-for aS = aS_values
-    % Format numbers
-    aS_str    = sprintf('%03.0f', aS);       
-    theta_str = sprintf('%03.0f', theta);    
-    phi_str   = sprintf('%03.0f', phi);      
-    N_str     = sprintf('%d', N);  
-
-    % Construct directory path
-    SaveDirectory = ['./Results/Data_3D/GradientDescent/Phi',theta_str, '/aS_', aS_str, ...
-                     '_theta_', theta_str, '_phi_', phi_str, '_N_', N_str];
-
-    % Call the plotting function
-    try
-        Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
-    catch ME
-        warning(ME.message)
-        continue;
-    end
-
-    % Pause to inspect plot before continuing
-    disp(['Plotted for aS = ', num2str(aS)])
-    pause(1.5)
-end
-
 %%
 SaveDirectory = './Results/Data_3D/PhaseDiagram/GradientDescent/';
 JobNumber     = 0;
 Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
 %%
-SaveDirectory = './Results/Data_3D/PhaseDiagram/ImagTimePropagation/aS_080_theta_000_phi_000_N_508333';
+SaveDirectory = './Results/Data_3D/PhaseDiagram/ImagTimePropagation/';
 JobNumber     = 0;
 Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
+%%
+SCATTERING_LENGTH_RANGE = [91.46 92.50 93.54 94.58 95.62 96.67 97.71];
+NUM_ATOMS_LIST          = [100000 304167 508333 712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 4591667 4795833 5000000];
+
+JobNumber = 0;
+
+for i = 1:length(SCATTERING_LENGTH_RANGE)
+    aS = SCATTERING_LENGTH_RANGE(i);
+    for j = 1:length(NUM_ATOMS_LIST)
+        N = NUM_ATOMS_LIST(j);
+        
+        SaveDirectory = sprintf('D:/Results/Data_3D/PhaseDiagram/aS_%.6e_theta_000_phi_000_N_%d', aS, N);
+        fprintf('Processing JobNumber %d: %s\n', JobNumber, SaveDirectory);
+
+        % Call the plotting function
+        try
+            Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
+        catch ME
+            warning(ME.message)
+            continue;
+        end
+    
+        % Pause to inspect plot before continuing
+        pause(1.5)
+    end
+end
 
 %% Visualize phase diagram
 load('./Results/Data_3D/GradientDescent/phase_diagram_matrix_theta_30.mat')
@@ -659,6 +643,68 @@ else
     disp('The state is not modulated');
 end
 
+%% Generate lists
+
+% Set display format
+format longG
+
+% Generate the lists
+list1 = linspace(1E5, 5E6, 25);
+list2 = linspace(80, 105, 25);
+
+% Convert to strings with no scientific notation
+str_list1 = compose('%.0f', list1);
+str_list2 = compose('%.2f', list2);
+
+% Join as space-separated strings
+row1 = strjoin(str_list1', ' ');
+row2 = strjoin(str_list2', ' ');
+
+% Display results
+disp(row1)
+disp(row2)
+
+%% Use space-separated floating-point/integer values
+SCATTERING_LENGTH_RANGE = "[80.00 81.04 82.08 83.12 84.17 85.21 86.25 87.29 88.33 89.38 90.42 91.46 92.50 93.54 94.58 95.62 96.67 97.71 98.75 99.79 100.83 101.88 102.92 103.96 105.00]";
+NUM_ATOMS_LIST          = "[100000 304167 508333 712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 4591667 4795833 5000000]";
+ 
+%% Phase diagram for untilted case
+N = [1E5, 3.04E5, 5.08E5, 7.125E5, 9.16E5, 1.12E6, 1.325E6, 1.529E6, ...
+     1.733E6, 1.9375E6, 2.141E6, 2.345E6, 2.55E6, 2.75E6, 2.95E6, ...
+     3.1625E6, 3.367E6, 3.57E6, 3.775E6, 3.979E6, 4.183E6, 4.3875E6, ...
+     4.591E6, 4.795E6, 5E6];
+
+as_UB = [88.33, 94.58, 95.62, 96.67, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, ...
+         97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71, 97.71];
+
+as_LB = [87.29, 93.54, 94.58, 95.62, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, ... 
+         96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67, 96.67];
+
+% Filter only rows with non-NaN data
+valid_idx = ~isnan(as_LB) & ~isnan(as_UB);
+N_valid   = N(valid_idx);
+LB_valid  = as_LB(valid_idx);
+UB_valid  = as_UB(valid_idx);
+
+% Create shaded area between LB and UB
+x_fill    = [N_valid, fliplr(N_valid)];
+y_fill    = [LB_valid, fliplr(UB_valid)];
+
+% Plot settings
+figure(1);
+set(gcf,'Position',[100 100 950 750])
+fill(x_fill, y_fill, [0.8 0.8 1], 'EdgeColor', 'none'); % Light blue shade
+% Axes settings
+set(gca, 'XScale', 'log');
+xlim([1E4, 1E7]);
+ylim([79, 106]);
+xlabel('Atom number', 'Interpreter', 'latex', 'FontSize', 16);
+ylabel('Scattering length $a_s$ ($a_0$)', 'Interpreter', 'latex', 'FontSize', 16);
+grid on;
+set(gca,'FontSize',16,'Box','On','Linewidth',2);
+
+%%
+
 function ModulationFlag = determineDensityModulation(psi, Params, Transf)
     
     % Axes scaling and coordinates in micrometers
@@ -704,59 +750,3 @@ function ModulationFlag = determineDensityModulation(psi, Params, Transf)
     end
 end
 
-%% Generate lists
-
-% Set display format
-format longG
-
-% Generate the lists
-list1 = linspace(1E5, 5E6, 25);
-list2 = linspace(80, 105, 25);
-
-% Convert to strings with no scientific notation
-str_list1 = compose('%.0f', list1);
-str_list2 = compose('%.2f', list2);
-
-% Join as space-separated strings
-row1 = strjoin(str_list1', ' ');
-row2 = strjoin(str_list2', ' ');
-
-% Display results
-disp(row1)
-disp(row2)
-
-%% Phase diagram for untilted case
-N = [1E5, 3.04E5, 5.08E5, 7.125E5, 9.16E5, 1.12E6, 1.325E6, 1.529E6, ...
-     1.733E6, 1.9375E6, 2.141E6, 2.345E6, 2.55E6, 2.75E6, 2.95E6, ...
-     3.1625E6, 3.367E6, 3.57E6, 3.775E6, 3.979E6, 4.183E6, 4.3875E6, ...
-     4.591E6, 4.795E6, 5E6];
-
-as_LB = [87.29, 93.54, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, ...
-         NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN];
-
-as_UB = [88.33, 94.58, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, ...
-         NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN];
-
-% Filter only rows with non-NaN data
-valid_idx = ~isnan(as_LB) & ~isnan(as_UB);
-N_valid   = N(valid_idx);
-LB_valid  = as_LB(valid_idx);
-UB_valid  = as_UB(valid_idx);
-
-% Create shaded area between LB and UB
-x_fill    = [N_valid, fliplr(N_valid)];
-y_fill    = [LB_valid, fliplr(UB_valid)];
-
-% Plot settings
-figure(1);
-set(gcf,'Position',[100 100 950 750])
-fill(x_fill, y_fill, [0.8 0.8 1], 'EdgeColor', 'none'); % Light blue shade
-% Axes settings
-set(gca, 'XScale', 'log');
-xlim([1E4, 1E7]);
-ylim([79, 106]);
-xlabel('Atom number', 'Interpreter', 'latex', 'FontSize', 16);
-ylabel('Scattering length $a_s$ ($a_0$)', 'Interpreter', 'latex', 'FontSize', 16);
-grid on;
-set(gca,'FontSize',16,'Box','On','Linewidth',2);
-